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1.  INTRODUCTION 


During  the  analysis  of  1 20-0)01  gun  firings  designed  to  look  at  tte  interior  ballistic  characteristics  of 
combustible  cartridge  cases  (Robbins,  Koszoru,  and  Minor  1986),  interior  ballistic  simulations  using 
XNOVAKTC  (XKTC)  (Gough  1980)  were  in  agreement  with  measured  pressure-time  curves  and 
pressure-difference  curves.  1BHVG2  (Anderson  and  Fickie  1987)  calculations  gave  a  maximum  breech 
pressure  that  was  42  MPa  higher  than  measured.  Parametric  studies  were  performed  using  XKTC  to 
attempt  to  attribute  this  disparity  to  the  various  processes  omitted  from  IBHVG2.  The  boattail  intrusion 
was  calculated  to  account  for  14  MPa,  with  effects  of  flamespreading  and  intergranular  stress  a^^counting 
for  3  MPa  each.  Subsequent  calculations  (Robbins  1986)  indicated  chambrage,  propellant  packaging,  wave 
dynamics,  and  multiphase  effects  (the  solid  propellant  velocity  lag  and  concomitant  fomialion  of  an  ullage 
region  between  the  projectile  base  and  the  propellant  bed)  as  contributors  to  the  differences  between  the 
lumped-parameter  and  two-phase  interior  ballistic  codes.  Further  study  demonstrated  that  the  influence 
of  chambrage  and  propellant  velocity  lag  could  be  represented  in  an  analytic  gradient  equation  (Robbins, 
Anderson,  and  Gough  1990).  This  report  extends  both  the  traditional  Lagrange  gradient  equation,  and  the 
chambrage  gradient  equation,  to  account  for  boattail  intrusion.  It  compares  the  analytic  pressure  gradient 
with  that  predicted  by  XKTC  and  assesses  the  extent  to  which  the  effects  of  boattail  intrusion  account  for 
the  differences  in  ballistic  predictions. 

The  family  of  NOVA  codes,  of  which  XKTC  is  the  latest  version,  has  been  used  with  uncompromised 
databases  to  model  gun  systems  successfully  (Robbins,  Koszoru,  and  Minor  1986;  Robbins  1983;  Robbins 
and  Horst  1984).  Since  XKTC  calculates  the  pressure  gradient  from  first  principles,  and  agrees  with  gun 
firings,  XKTC  simulations  are  assumed  correct.  Accordingly,  all  the  lumped-parameter  computer  runs, 
with  the  different  boattail  gradients,  are  compared  with  equivalent  XKTC  computer  runs. 

1.1  Models.  The  analysis  of  the  chambrage  gradient  equation  can  be  traced  back  as  far  as  Vinti 
(1942).  The  analyses  of  the  chambrage  gradient  equation,  the  original  analyses  of  the  propellant  velocity 
lag  gradient  equation,  and  the  combination  of  the  two  were  done  by  Gough  (Robbins,  Anderson,  and 
Gough  1990),  who  is  also  responsible  for  the  initial  development  of  the  boattail  gradient  equation  (Gough 
in  preparation).  Similar  analysis  for  the  chambrage  gradient  for  a  different  purpose  has  also  been 
performed  (Morrison  and  Wren  1989).  Reasonable  assumptions  for  the  accommodation  of  a  boattail 
intrusion  were  made,  and  this  document  represents  the  complete  analysis. 
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1.2  Influence  of  Chambrage  With  a  Boattail  Intrusion.  For  the  chambrage  gradiertt  equation  with 
boatiail  e'Tects,  the  propellant  is  assumed  to  be  uniformly  distributed  between  the  breech  and  the  base  of 
the  projectile.  The  variation  in  tube  area  is  assumed  to  be  confined  to  the  chamber,  while  the  area  of  the 
bore  is  constam.  The  boattail  is  assumed  to  be  a  right  circular  cylinder.  This  keeps  die  solution  consistent 
with  Uk;  simplicity  of  the  lumped-parameter  code,  reduces  computation  time  (cuts  down  on  integration) 
and  greatly  simplifies  the  derivation. 


The  continuity  and  momentum  equations  for  the  unsteady  flow  of  a  homogeneous,  inviscid  substance 
through  a  tube  with  variable  area  are 


a(Ap)  ^  aCpAu)  .  ^ 

sr 


(1) 


^  3u  ^  du  dP 

p-,—  +  pUi—  +  *  0  ,  (2) 

dt  dz  dz 

where 

A  =  cross-sectional  area 
P  =  pressure 
p  =  density 
u  =  velocity 
t  =  time 
z  =  distance. 

The  system  to  be  modeled  is  shown  in  Figure  1.  With  reference  to  the  figure, 

za  =  distance  from  breech  to  afit  end  of  projectile 

zp  =  distance  from  breech  to  base  of  projectile 

Ay^  =  cross-sectional  area  of  the  boattail 

Aqa  =  cross-sectional  area  of  projectile  base  (excluding  boattail) 

Ag  =  cross-sectional  area  of  bore 

^BA 

A  =  cross-sectional  area  of  chamber 

Aj  =  cross-sectional  area  of  the  tube  in  the  area  of  the  boattail 
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A  =^  +  AA 


also 


Vp  =s  projectile  velocity 
V(z,t)  =  volume  up  to  position  z,  at  time  t 
A(z,t)  =  area  at  position  z,  at  time  t 
V(zp)  =  volume  up  to  base  of  projectile. 


The  analysis,  a  detailed  version  of  which  is  shown  in  Appendix  A,  considers  two  regions  separately,  the 
region  from  the  breech  to  the  aft  end  of  the  projectile  and  the  region  from  the  aft  end  of  the  projectile  to 
base  of  the  projectile.  The  analysis  starts  with  the  Lagrange  assumption 


ap 

■ST 


-  0  . 


For 


0  :s  z  za  , 


the  velocity,  from  the  continuity  equation  (1),  is 


ABVpV(z.t) 
V(zp)A(z,t)  ’ 


(4) 
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and  the  velocity  for 


za  <  z  ^  zp 


is 


ABVpV(z,t) 

V(zp)A(z.t) 


AaVp 

A(z.t) 


(5) 


For  the  portion  of  space  behind  the  projectile,  0  ^  z  <  za.  calculations  of  the  pressure  distribution  from 
the  breech  to  the  aft  end  of  the  projectile  can  be  performed  by  taking: 

P(z,t)  =  pressure  at  z  and  t 

P(za,t)  =  P(za)  =  pressure  at  aft  end  of  projeaile  at  time  t 

C  =  total  charge  mass 

PgB  =  pressure  at  the  breech 

Pg  =  pressure  at  the  base  of  the  projectile 

Pres  =  resistive  pressure  to  the  motion  of  the  projectile 

Mp  =  mass  of  the  projectile, 

projectile  acceleration. 


V 


p 


P(za)Ay^  +  PgAgA  "  ArPres 


(6) 


then  using  (4)  and  the  momentum  equation  (2),  P(za)  can  be  written  as 


P(za)  -  zagfOPgR  +  za,(i)  +  zajCOPg , 


(7) 


where 


zao(t)  - - 

j  ^  CAgAy^Qjfza) 

v2(zp)Mp 


4 


CA(jVpQj(/a)  ^  CAgPjjg^QjCza)  CAgVpQ2(za) 

V\zp)  v2(zp)M_  2v\zp) 

za,(t)  =  - E - - - 

r  CAgAyyQj(za) 

1  +  - 

v2(zp)Mp 


za2(t)  = 


CAgAgys^QjCza) 

v2(zp)Mp 

CAgAAQi(za) 
+  — 
v2(zp)Mp 


and 


Q,(za) 


J 


A(za  ■) 


Q2(za) 


V^(za~) 
A  ^(za  ~) 


where 


za  ■  =  Jimj.^(,za  -  e  . 


The  pressure  distribution  P(z)  is 


P(z) 


HR 


a3(t)  +  a5(t)Pgg  +  aj(t)Pg  jQ]  (z)  +  b(t)Q2(z), 


where 


a3(t) 


2  2 
CAg^Vp 

V\zp) 


CAgA^za,(t)  ^  CAgPggj 
v2(zp)Mp  v2(zp)Mp  ’ 


(8) 
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v2(zp)Mp"  v2(zp)Mp' 


ajd)  =  - 


v2(zp)Mp 


b(t)  =  - 


2V3(zp)’ 


and 


Q2(^-) 


vV-) 

A  Vo 


Now,  focusing  on  the  portion  from  the  aft  end  of  the  projectile  to  the  base  of  the  projectile... 


za  ■<  z  i  zp  , 


P(z)  =  zao(t)PBR  +  zai(t)  +  zajOjPg  +  fk  +  ajfOQjfza  *.t) 

+  a4(t)PBQ,(za  *,t)  +  a5(t)PBRQ,(za  *,t)  +  CjfOQjCza  *,t) 

+  C4(t)PBQ3(za  \t)  +  C5(t)PBRQ3(za  *.l)  bfOQjfz.t) 

+  h,Q4(z,t)  +  j,Q5(z,t)  +  k, . 


(9) 
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where 


CVpAaAb 

V2(zp) 


Jl  = 


cA^v; 

2V(zp) 


za  *  -  Um^_^.a  +  e  , 


CVpA^V^za*)  _  CVp^A^A|,V(/-a-)  ^  CA>^^ 

2V\zp)A^(za*)  v2(zp)A2(za*)  2V(zp)A  ^(za  *)  ’ 


Qa(z.t)  = 


A^(z,t) 


QsCai)  = 


aV-O 


Q^(za^t)  =  jZ^dz , 
A(z.l) 


Q3(za*,l)  =  J* 


dz 


A(z.l) 


za 
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C3(t) 


CA^zai(t) 

V(zp)Mp 


CAy^AaPRES 

V(zp)Mp 


CVpA^A, 

V2(zp) 


and 


C4(t) 


^  CA^zajd) 
V(zp)Mp  V(zp)Mp  ’ 


Cjd) 


CA^zaod) 
V(zp)Mp  ■ 


fk  =  jump  in  pressure  across  the  boattail 

An  analysis  by  Kooker  (April  1991)  indicates  that  the  pressure  drop  across  the  boattail  (which  in  a 
one-dimensional  analysis  is  equivalent  to  determining  the  pressure  drop  across  a  moving  discontinuity  in 
area,  Figure  2)  is  given  by: 


Figure  2.  Presentation  of  a  moving  area  in  tube. 


Mass  balance: 


Pi(“i  -  *  P2(“2  -  'Vp)A2  =  rti 


(10) 
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Momentum  balance: 


P2U2(U2  -  Vp)A2  -  PiUjCu,  -  Vp)A,  *  PjAj  -  P2A2  >  Pmean(A2 


(11) 


=  -(P2 


-Pi) 


A  j  A2 


since 

=  *^2) 


Thus, 


where 


(P2  -  Pi) 


2p(u,  -  Vp)A,(U2  -  u,) 
A  j  +  A2 


A,  =  A(za‘), 


(12) 


A2  =  A(/,a*). 


Pi  =  P2  = 


C 

V(zp)  ’ 


u,  =  u(/.a"). 


U2  =  u(/.a*). 


or 


-2 

r  c  1 

,Vp^(7..-)A, 

ABV(za)  ^ 

2 

[V(zp)  J 

*  A(za  *) 

V(zp)A(za-) 

(A(za‘)  +  A(za*)) 


(13) 


9 


Wc  arc  now  at  a  point  whcrc  the  projectile  base  and  breech  prcssurc  can  be  determined.  For  use  in 
lumped-parameter  interior  ballistic  models,  the  gradient  equation  is  usually  cast  in  terms  of  the  mean 
prcssure  (P„): 


/.p 

jA(z.t)P(z.t)3/ 

0 

zp 

jA(z.t)dz 


Substituting  the  pressure  distribution  and  integrating. 


P^  =  b,(zp.t)  +  bjCzp.OPg  +  bjCzp.tjPgg  . 


bi(/,p,t)  = 


a3(t)Q7(za)  +  a3(t)Q9(zp)  +  b(t)Q6(zp) 
_ 


za,(t)(V(zp)  -  V(za))  +  nc(V(zp)-V(za)) 
_ 

C3(1)(;  (z->)  +  h,Q,(zp)  +  j,Q3(zp) 

V(zp) 


k,(V(zp)  -  V(za)) 
V(^P) 


b2(zp.t)  = 


a4(t)Q7(za)  +  a4(t)Q9(zp)  ^  za^(t)(V(zp)  -  V(za)) 


V(zp) 


C4(t)Q8(zp) 

V(zp) 
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b3(zp.t) 


V(za)  +  35(007(23)  +  35(009(2?) 
2ao(tKV(2p)  -  V(za))  +  CjOgC^P) 

*  V(^) 


V2(z,t) 

A(z.i) 


•dz , 


and 


07(23)  =  r  A(z.O  . 


Q,<zp)  .  jAfe.1) 

za  *  za  * 


09(2?)  = 


V(z.O 

A(z,0 


d  /d  v. . 


P(zp)  can  be  dcicrmincd  from  (9)  and  can  be  wriitcn  as 


P(2P)  =  Pb 


l3(2P,0p  ^  l,(zp,0 

1;^  ^  ’ 


where 


l,(zp.O  =  za,(0  +  fie  +  a3(00i(zp)  +  03(003(2?) 

+  b(t)02(zp)  +  h,Q4(zp)  +  j,05(2p)  +  k, , 


(16) 


II 


l2(/.p.i)  »  1  -  zajd)  -  a4(t)Q,(zp)  -  C4(t)Q3(zp) . 
I3(zp.i)  »  zao(l)  +  35(1)0, (zp)  +  C5(t)Q3(zp) , 


Qi(zp) 


f  V(z.t) 
J.A(z.t) 


dz , 


QjCzp) 


V^(2p) 

aV-p)  ’ 


zp 

Q3(ZP)  =  J- 


dz 


A(z.t) 


Q4(>'-P) 


V(zp) 

aVp)’ 


and 


QjCzp)  = - 

A  ^(zp) 


Therefore,  (15)  and  (16)  arc  solved  simultaneously: 


^  b,(zp.t)  ^  l,(zp,i) 

b3(zp,i)  b3(zp,i)  l3(zp,t) 

b2(zp,t)  l2(zp,t) 

+  ___ 

b3(zp.t)  l3(zp,t) 


(17) 
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and 


(18) 


l2(zp.t)p  li(zp.t) 

l3(zp.t)  ®  l3(zp.t) 

The  energy  in  the  fluid  is  represented  by 

dE  »  J-u^dm  ,  (19) 

2 

dm  «  pAdz  , 

and  using  (4)  and  (S), 

zp 

E  =  JdE  =  -  b(t)Q6(zp)  -  h,Q,(zp)  -  j,Q3(zp)  .  (20) 

0 

2.  CALCULATIONS 

Interior  ballistic  simulations  were  performed  with  IBRGAB  (Appendix  B — an  interior  ballistic  code 
into  which  the  boattail  gradient  had  been  incorporated)  to  investigate  the  influence  on  the  interior  ballistic 
trajectory  of  a  flat-based  projectile  with  a  boattail  (with  and  without  chambrage).  The  same  simulations 
were  done  with  XKTC  using  XKTC  databases  as  consistent  with  those  of  IBRGAB  as  the  physical  scope 
of  XKTC  allows. 

The  calculations  performed  with  IBRGAB  involved  databases  with  evenly  distributed,  seven -perforated 
propellant  having  an  initial  porosity  of  0.4579  and  zero  barrel  resistance,  and  assumed  the  propellant 
ignited  at  an  initial  instant.  All  calculations  were  performed  for  a  flat-based  projectile  with  and  without 
a  boattail  and  no  heat  loss. 

The  parameters  used  in  the  computer  codes  were: 

Bore  diameter  127  mm 

Volume  .01  m^ 

Travel  4.572  m 

Propellant  mass  9.0  kg 

Projectile  mass  2.25,  9.0,  36.0. 
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Propellant  characteristics  were: 


Impetus 

Covolume 

Gamma 

Rame  temperature 
Molecular  weight 
Density 
Burning  rate 


.1136  X  lO’j/kg 
.976  X  10"^  m^Ag 
1.23 

3,143  K 

23.0  kg/kg  mole 
.16605  X  lO^kg/m^ 

.001 10519P‘  °  m/s  (P  is  in  MPa). 


The  maximum  breech  pressure  studied  in  IBRGAB  was  500  MPa.  It  was  obtained  by  varying  the 
outer  diameter  of  the  propellant  grain,  while  holding  the  grain  length  constant  at  .03175m.  Calculations 
were  performed  for  a  combination  of  four  chambers,  each  with  four  boattaii  values: 

1)  0%  chambrage  with  0.  5.  10.  15%  boattaii 

2)  5%  chambrage  with  0.  5.  10.  15%  boattaii 

3)  10%  chambrage  with  0,  5,  10,  15%  boattaii 

4)  15%  chambrage  with  0,  5,  10,  15%  boattaii. 

The  chamber  with  0%  chambrage  and  0%  boattaii  is  a  straight  tube  with  the  chamber  diameter  equal  to 
the  bore  diameter.  The  5,  10,  and  15%  chambrage  chambers  were  obtained  by  adding  5, 10,  and  15%  of 
the  bore  diameter,  respectively,  to  the  diameter  of  the  breech  for  the  0%  chambrage  chamber.  With  the 
increase  in  chambrage,  there  is  a  decrease  in  chamber  length  to  keep  the  volume  constant.  The  breech 
diameter  increases  from  .127  m  to  .146  m  (0  to  15%),  while  the  chamber  length  decreased  to  .78430  m 
from  .90782  m.  The  boattaii  length  is  held  constant  at  .508  m  and  the  boattaii  percentages  sim{4y 
represent  that  fraction  of  the  constant  chamber  volume.  The  different  chambrage  and  boattaii 
combinations  are  arranged  in  tables  according  to  the  ratio  of  the  charge  mass  to  the  projectile  mass  (c/m). 
In  the  tables,  the  nomenclature,  bt.  refers  to  the  boattaii  aixi  ch  refers  to  chambrage.  Each  maximum 
breech  pressure  and  muzzle  velocity  drop  is  clearly  seen  in  both  XKTC  and  IBRGAB  portions  of  the 
tables.  Table  1  gives  tire  baseline  maximum  breech  pressures  and  muzzle  velocities  for  each  c/m  in 
XKTC  and  IBRGAB.  In  IBRGAB,  maximum  breech  pressure  is  forced  to  5{X)  MPa  by  varying  the  grain 
diameter.  An  equivalent  XKTC  database  was  then  made  to  yield  the  values  given.  Tables  2  through  4 
contain  the  explicit  results  of  the  different  combinations  of  boattaii  and  chambrage,  emphasizing  the  effects 
of  the  boattaii  for  c/m  of  .25, 1.0,  and  4.0.  For  example,  in  any  table,  to  go  from  the  baseline  calculation 
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to  a  10%  boattail  and  a  5%  chambrage.  a  boattail  is  added  which  is  10%  of  the  chamber  volume,  and  the 
chamber  description  is  altered  by  adding  S%  of  the  bore  diameter  to  the  diameter  of  the  breech.  Arranged 
in  the  same  form  are  taUes  S  through  6  depicting  the  differences  from  the  effects  of  the  boattail.  The 
differences  are  obtained  by  subtracting  either  the  S,  10.  or  1S%  boattail  result  from  the  baseline  result. 
These  values  allow  each  table  to  clearly  illustrate  that  a  larger  boattail  combined  with  a  larger  chambrage 
give  a  greater  drop  in  maximum  breech  pressure  and  to  show  the  maximum  breech  pressure  drop  increases 
while  the  cAn  ratio  increases. 

3.  RESULTS 

3.1  Baselines. 


Table  1.  Baselines  (0  Boattail.  0  Chambrage) 


C/M 

XKTC 

IBRGAB 

— 

**MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(nVs) 

**MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

O.D. 

(m) 

.25 

502.6 

898.8 

500.0 

899.5 

.016599 

1.0 

540.3 

1,641.3 

499.9 

1,596.8 

.009879 

4.0 

512.8 

2,538.4 

499.9 

2,342.6 

.007237 

3.2  IB  Calculations  With  Boattail  and  Chambrage. 


Table  2.  IB  Calculation  With  Boattail  and  Chambrage,  C/M  of  .25 


OUTPUT 

XKTC 

IBRGAB 

%bt 

%ch 

f*MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

**MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

00 

00 

502.6 

898.8 

500.0 

899.5 

05 

00 

501.0 

897.9 

497.9 

898.9 

10 

00 

496.1 

896.8 

495.9 

898.3 

Tabic  2.  IB  Calculation  With  Boatlail  and  Chambrage,  C/M  of  .25  (continued) 


OUTPUT  XKTC  IBRGAB 


%bt 

%ch 

^’max 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

**MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

15 

00 

497.8 

896.1 

494.0 

897.8 

(X) 

05 

500.6 

898.2 

497.8 

898.8 

05 

05 

499  2 

897.6 

495.5 

898.1 

10 

05 

496.0 

895.7 

493.4 

897.4 

15 

05 

498.8 

897.1 

491.4 

896.8 

00 

10 

498.8 

897.8 

495.8 

898.1 

05 

10 

497.2 

897.2 

495.2 

898.1 

10 

10 

494.8 

895.9 

491.0 

896.7 

15 

10 

489.4 

893.2 

488.8 

896.0 

(X) 

15 

497.0 

897.3 

494.1 

897.6 

05 

15 

495.7 

896.8 

491.3 

896.7 

10 

15 

493.3 

896.0 

488.7 

896.0 

15 

15 

489.0 

893.6 

486.3 

895.2 

Table  3.  IB  Calculation  With  Boattail  and  Chambrage,  C/M  of  1.0 

OUTPUT 

XKTC 

IBRGAB 

f^MAX 

MUZZLE 

**MAX 

MUZZLE 

%bt 

%ch 

BREECH 

VELOCITY 

BREECH 

VELOCITY 

(MPa) 

(m/s) 

(MPa) 

(m/s) 

00 

00 

540.3 

1,641,3 

499.9 

1,596.8 

05 

00 

530.2 

1,635.2 

493.4 

1,593.1 

10 

00 

516.3 

1 ,626.9 

487.1 

1,589.5 

15 

00 

509.5 

1,623.7 

481.3 

1,586.0 

00 

05 

531.7 

1,638.3 

492.8 

1,592.0 
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Table  3.  IB  Calculation  With  Boattail  and  Chambrage,  C/M  of  1.0  (continued) 


OUTPUT 

XKTC 

IBRGAB 

%bt 

%ch 

^MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

f*MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

05 

05 

523.0 

1.632.7 

485.7 

1,587.8 

10 

05 

510.8 

1,625.2 

479.0 

1,583.6 

15 

05 

507.0 

1,623.2 

472.7 

1,579.8 

00 

10 

523.5 

1,635.3 

486.4 

1,587.8 

05 

10 

515.0 

1,629.9 

480.3 

1,584.4 

10 

10 

505.0 

1.623.1 

471.3 

1,578.7 

15 

10 

491.0 

1,624.6 

464.7 

1.574.2 

00 

15 

515.8 

1,632.3 

480.9 

1.584.3 

05 

15 

507.9 

1,626.9 

472.3 

1,579.1 

10 

15 

498.3 

1,620.5 

464.4 

1.574.0 

15 

15 

486.6 

1,612.1 

457.1 

1.569.0 

Table  4.  IB  Calculation  With  Boattail  and  Chambrage,  C/M  of  4.0 


OUTPUT 

XKTC 

IBRGAB 

%bt 

%ch 

^MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

f^MAX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

00 

00 

512.8 

2.538.4 

499.9 

2,342.6 

05 

00 

510.1 

2,523.1 

485.6 

2,328.5 

10 

00 

505.9 

2,499.9 

472.3 

2,314.7 

15 

00 

507.5 

2.486.7 

460.0 

2,301.1 

00 

05 

503.0 

2,502.5 

484.9 

2,323.8 

05 

05 

496.3 

2,484.9 

469.6 

2,. 307. 4 

10 

05 

489.6 

2,461.1 

455.5 

2,290.9 
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Tabic  4.  IB  Calculation  With  Boaltail  and  Chambragc,  C/M  of  4  0  (continued) 


OUTPUT 

XKTC 

IBRGAB 

%bt 

%ch 

PmaX 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

Pmax 

BREECH 

(MPa) 

MUZZLE 

VELOCITY 

(m/s) 

15 

05 

490. 1 

2,448.9 

442.5 

2,276.9 

(X) 

10 

494.0 

2,469.9 

471.7 

2,306.7 

05 

10 

484.2 

2,449.0 

456.6 

2,290.1 

10 

10 

476.1 

2,424.0 

439.9 

2,270,0 

D 

10 

464.3 

2,391.7 

426.2 

2,252.3 

(X) 

15 

485.1 

2,440.3 

460.3 

2,291.8 

05 

15 

472.8 

2,415.9 

442.2 

2,270.6 

10 

15 

463.6 

2,389.1 

425.8 

2,250.3 

15 

15 

452.1 

2,356.9 

411.1 

2,2.30.8 

3.3.  Droo  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail, 


Table  5.  Drop  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail,  C/M  of  .25 


-  .  . 

OUTPUT 

XKTC 

%bi 

%ch 

APmaX  BRCH 
(MPa) 

AMUZZLE  VEL 
(m/s) 

(X) 

(X) 

0.0 

0.0 

05 

(X) 

1.6 

0.9 

10 

00 

6.5 

2.0 

15 

(X) 

4.8 

2.7 

00 

05 

2.0 

0.6 

05 

05 

3.4 

1.2 

10 

05 

6.6 

3.1 

15 

05 

3.8 

1.7 

(X) 

10 

3.8 

1.0 

IBRGAB 


APmaX  BRCH 
(MPa) 


AMUZZLE  VEL 
(m/s) 


0.0 


0.0 


2.1 


0.6 


4.1 


1.2 


6.0 


1.7 


2.2 


0.6 


4.5 


1.4 


6.6 


2.1 


8.6 


3.3 


4.2 


1.4 


Tabic  5.  Drop  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail,  C/M  of  .25  (continued) 


XKTC 


IBRGAB 


AMUZZLE  VEL. 
(m/s) 

aPmax  bRCh 
(MPa) 

1.6 

4.8 

2.9 

9.0 

5.6 

11.2 

1.5 

5.9 

2.0 

8.7 

2.8 

11.3 

5.2 

13.7 

Table  6.  Drop  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail,  C/M  of  1.0 


OUTPUT 

XKTC 

IBF 

%bt 

%ch 

^Pmax  bRCH 
(MPa) 

AMUZZLE  VEL. 
(m/s) 

APmaX  BRCH 
(MPa) 

00 

00 

0.0 

0.0 

0.0 

05 

00 

10.1 

6.1 

6.5 

10 

00 

24.0 

14.4 

12.8 

15 

00 

30.8 

17.6 

18.6 

00 

05 

8.6 

3.0 

7.1 

05 

05 

17.3 

8.6 

14.2 

10 

05 

29.5 

16.1 

20.9 

15 

05 

33.3 

18.1 

27.2 

00 

10 

16.8 

6.0 

13.5 

05 

10 

25.3 

11.4 

19.6 

10 

10 

35.3 

18.2 

28.6 

15 

10 

49.3 

16.7 

35.2 

AMUZZLE  VEL. 
(m/s) 
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Tabic  6.  Drop  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail,  C/M  of  1.0  (continued) 


OUTPUT 

XKTC 

IBRGAB 

%bi 

%ch 

APmax  BRCH 
(MPa) 

AMUZZLE  VEL. 
(m/s) 

aPmax  brch 

(MPa) 

AMUZZLE  VEL 
(m/s) 

(X) 

15 

24.5 

9.0 

19.0 

12.5 

05 

15 

32.4 

14.4 

27.6 

17.7 

10 

15 

42.0 

20.8 

35.5 

22.8 

15 

15 

53.7 

29.2 

42.8 

27.8 

Table  7.  Drop  in  Pressure  and  Muzzle  Velocity  Due  to  Boattail,  C/M  of  4.0 


OUTPUT 

XKTC 

IBRGAB 

%bl 

%ch 

APn^ax  brch 

(MPa) 

AMUZZLE  VEL. 
(m/s) 

aPmax  brch 

(MPa) 

AMUZZLE  VEL. 
(m/s) 

00 

00 

0.0 

0.0 

0.0 

0.0 

05 

00 

2.7 

15.3 

14.3 

14.1 

10 

00 

6.9 

38.5 

27.6 

27.9 

15 

00 

5.3 

51.7 

39.9 

41.5 

00 

05 

9.8 

35.9 

15.0 

18.8 

05 

05 

16.5 

53.5 

30.3 

35.2 

10 

05 

23.2 

77.3 

44.4 

51.7 

15 

05 

22.8 

89.5 

57.4 

65.7 

00 

10 

18.8 

68.5 

28.2 

35.9 

05 

10 

28.6 

89.4 

43.3 

52.5 

10 

10 

36.7 

114.4 

60.0 

72.6 

15 

10 

48.5 

146.7 

73.7 

90.3 

00 

15 

27.7 

98.1 

39.6 

50.8 

05 

15 

40.0 

122.5 

57.7 

72.0 

10 

15 

49.2 

149.3 

74.1 

92.3 

15 

15 

60.7 

181.5 

88.8 

111.8 

Il  was  expected  that  the  introduction  of  a  boattail  without  any  volu.  ■<€  change  would  lead  to  a  lowering 
of  the  maximum  breech  pressure.  This  was  expected  because  the  aft  end  of  the  projectile  is  subject  to  a 
higher  pressure  than  the  base  of  the  projectile  and  therefore  is  accelerated  faster.  More  chamber  volume 
is  then  opened  up.  leading  to  a  decrease  in  maximum  breech  pressure  (P^ax)- 

In  the  baseline  analysis  (Table  1, 0  boattail,  0  chambrage),  the  calculated  maximum  breech  pressures 
are  within  10%  of  the  maximum  breech  pressures  given  by  XKTC  for  a  comparable  database,  and  the 
velocities  arc  within  8%. 

With  a  c/m  of  .25,  IBRGAB  captures  the  boattail  effects  as  measured  by  the  change  in  maximum 
breech  pressure  and  velocity  in  a  qualitative,  if  not  quantitative  manner.  When  compared  to  XKTC,  the 
calculations  with  15%  boattail  and  chambrage  of  0  and  5%  exhibit  a  rise  in  maximum  breech  pressure 
above  that  with  10%  boattail.  This  was  unexpected. 

Again,  for  a  c/m  of  1.0,  the  agreement  in  the  change  in  maximum  breech  pressure  is  good.  Note  that 
the  change  in  maximum  breech  pressure  in  IBRGAB  is  slightly  smaller  than  the  change  in  maximum 
breech  pressure  in  XKTC,  yet  both  are  rather  uniform. 

La.sily,  at  c/m  of  4.0,  the  pressure  drop  in  XKTC  for  small  boaiiails  is  unexpectedly  small  (less  than 
for  c/m  =  1).  It  is  al.so  .smaller  than  that  of  the  analytic  gradient  with  boattail.  Similar  effects  ol  the 
reversal  of  the  change  in  maximum  breech  pressures  were  observerd  for  a  c/m  of  .25  in  XKTC. 

4.  CONCLUSIONS 

A  boattail  has  been  incorporated  into  the  gradient  equation  of  a  lumpcd-paramcicr  interior  ballistics 
model  with  reasonable  effccLs. 

The  gradient  equation,  with  boattail  addition,  captures  qualitatively,  if  not  quantitatively,  the  effects 
for  normal  c/m’s  (.25  and  1.0). 

For  c/m  of  4.0,  or  large  c/m’s,  some  questions  remain  about  the  physics  of  the  gradient  model. 
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APPENDIX  A: 

THE  DERIVATION  OF  A  GRADIENT  EQUATION 
WITH  AREA  CHANGE  IN  A  TUBE  AND  A  BOATTAIL 
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The  continuity  and  momentum  equations  for  the  unsteady  flow  of 
a  homogeneous,  inviscid  substance  through  a  tube  with  variable  area 
are 


d(A(z,  t)  p)  ^  d(pA(z,  t)  u(z,  t) )  ^  Q 
dt  dz 


du{z,  t) 


du(z,  t) 


where  A  =  cross  sectional  area 
P  =  pressure 
p  =  density 
u  =  velocity 
t  =  time 
z  =  distance. 

If  we  take  the  boattail  to  be  a  right  circular  cylinder  (Figure 
A-1)  , 


then  Ag  =  Ag;;^  +  Ap^,  where 

A^  =  cross-sectional  area  of  the  boattail, 

•^BA  ~  Area  of  the  base  of  the  projectile  exposed  to  fluid, 

Ag  =  Bore  area; 

and  A  =  A;^  +  A^,  where 

A  =  external  area  of  the  tube, 

A^  =  internal  area  of  the  tube (account  for  boattail). 

We  start  by  examining  the  area  from  the  breech  to  the  aft  end  of 
the  projectile. 

0  z  <  za  , 


where  A  is  the  area  associated  with  the  chamber  and  tube  wall. 
Performing  the  indicated  differentiation  of  (1)  and  making  the 
Lagrange  assumption. 


(1)  becomes 


dA(z,  t)  u(z,  t)  _  A{z,  t) 

dz  p  dt 


dA{z,  t) 
Bt 


The  density  can  be  written  as 

P  {of  chef  luid) 


V{zp) 


(3) 


(4) 


The  fluid  is  considered  to  be  composed  of  solid  and  gaseous 
propellai.t 

C  =  total  mass  of  propellant  (fluid), 
or 

C  -  initial  mass  of  propellant  (solid), 
and 

V(zp)  =  chamber  volume  up  to  the  base  of  the  projectile. 

Differentiating  the  density  gives 

Bp  _  p  Bv(zp) 

Bt  V(zp)  Bt 


but 


BV{zp) 

Bt 


where  Vp=velocity 
With  the  boundary 


of  the  projectile, 
condition 

L/(0)  =  0 


and 


BA{z,  t) 
Bt 


0 


for 


0  i  z  <  za  , 


28 


(3)  becomes 


u(z) 


A^V^V{z.  t) 
A{z,  t)  V(zp) 


(5) 


...from  the  aft  end  of  the  boattail  to  the  projectile  base, 

za  <  z  i.  zp  . 


Once  again  performing  the  indicated  differentiation  of  (1)  and 
making  the  Lagrange  assumption. 


d(A(z,  t)  u(z,  t)  )  _  A(z,  t)AgVjp  dA(z,  t) 
dz  V(zp)  dt 


Since  the  boattail  is  a  right  circular  cylinder 

0A(z,  t)  _ 

at  ' 


where  A  =  accounts  for  boattail  area; 
for 

za  <  z  i  zp 


and 

Z  Z 

J  d(A(z,  t)  u(z,  C)  )  -  j"  A(z,  t:)dz  , 

za*  za' 


where 


za*  =  +  e  . 


A{z,t)u{z,t)  -  A{za*)  u(za*)  =  [v(z,  t)  -V{za*,t)]  , 

V{zp) 

and  with  the  boundary  condition 

uizp)  =  Vj,  , 
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{V[zp)A{zp)  -  Ag(V(zp)  -  V(za*)) 
V(zp)  A{za*) 


u(za*)  . 


Since  Ag  =  +  Ag;^  and  A(2p)=Ab;^, 


u{za*) 


i 


A^V(za*)  -Aj^V(zp)\ 
v(zp)A(za*)  j 


Therefore, 


A(z.t)u(z.t)  =  {V(z.  t)  -  V{za*))  +  vj- 

V{zp) 


V{zp) 


and 


u(z,  t) 


V^^Viz.C)  _ 
V(zp)A(z,  t)  A{z,  t) 


(6) 


To  obtain  the  pressure  distribution,  we  once  again  examine  the  area 
from  the  breech  to  the  aft  end  of  the  boattail. 

0  i  z  <  za  . 


Differentiate  (5)  with  respect  to  time. 


du  _ 

t) 

dt 

v\zp)A{z,  t) 

V(zp)A{z.t) 

AgVj,V(Z, 

t)  ldv{,zp)\ 

AgVpV(z,t)  idA{z,t)\ 

( zp)  A  ( Z 

,t)l  at  / 

v(zp)A^{z,t)\  at  / 

with 
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dVjz,  t) 
dt 


dV(zp) 

dt 


Ba{z,  t) 
dt 


=  0  , 


=  0  , 


then 


du{z,  t)  ^  A^VpViz,  t)  _  A^V^Viz,  t) 
dt  V(zp)A(.z.t)  V^[zp)A{z,t) 


Differentiate  (5)  with  respect  to  distance. 


du(z,  t) 

A^V^V(z,t) 

ldA(z,  t) 

dz 

V{zp)A{z,  t) 

V(zp)A^lz,  t) 

\  dz 

with 


dV{z,  t) 
dz 


A(z,  t)  , 


then 


du(z,  t) 

A^V^V{z,t) 

idAiz,  t)  \ 

dz 

V{zp) 

V{zp)A^{z,  t) 

dz  j 

Substitute  (4),  (5),  (7),  and  (8)  into  equation  (2),  the 

equation,  and  we  get  the  equation: 


(8) 

momentum 
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Or 


c 

'  AgV^Viz.t) 

V(zp) 

^  V{zp)A{z,  t) 

Aiv^V(z,t)  ] 
V^{zp)A{z,  t)  j 


AgV^V(z.C)  ( 

^bV,  \ 

V(zp) 

V(zp)A{z.  t)  ( 

v{zp)  ) 

C 

AgV^Viz.t)  / 

AgV^V(z.t) 

ldA{z,  t) 

v(zp) 

V(zp)A(z,  t)  ( 

v(zp)A^(z.  t) 

\  dz 

dP  ^ 

-CA^V^V(z.t) 

CAiV^VHz.t)  iaA(z,t)\ 

dz 

V^{zp)A{z,  t) 

{zp)A^  (z.  t)\  3z  / 

and 

hp  =  f  C)  a.  ,  r  alv^vuz.t)  ,aa(z,t)u_. 

{  {  (zp)A(z,  t)  {  izp}A^  {z.  t)\  dz  j 


With  the  definition  P(0)  =  Pgg 


P(z)  =  P, 


Br 


^B^p  f(  V(Z,  t) 
V^(zp)  {xAiz,  t) 


+  rv^iz.t)  afl(z,  t)  g„ 

vUzp)  I  A^(z,  t)  dz 


(9) 


Using  Integration  by  Parts  for  the  second  part  of  equation  (9), 

Judv  =  uv  -  jvdu 


with  the  following  substitutions: 

dv  = 

A^iz,  t) 
u  =  V^(z,  t) 

V  - - - - - 

2A{.z,  t)* 

du  =  2V{z,  t)  dv  =  2V(z,  t)A(z,  t)dz 
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C  v^(z,  t)  dA{z,  t)  (z,  t)  ^  r  2V(z,  t)A(z,  t) 

JAUz,t)  dz  ~  2A^{z.t)  J  2AHz,t) 


r  V^iz,  t)  dA(z,  t)  _  -V^  iz,  t)  }  V(z,  t)  a, 
lAUz,t)  dz  '^^-2AUz,t)  {Aiz.t)"^^- 


By  substituting  from  equation  (10)  into  equation  (9)  and  factoring, 
the  result  becomes 


P(z) 


=  P, 


Br 


CA^Vy  ^  ^bVq  ‘fV(z.t) 
V^(zp)  {zp)  \{  A(z.  t) 


CA^Vp  VUz,t) 
2V^(zp)  A^(z.  t) 


dz  . 


(11) 


The  acceleration  of  the  projectile,  dVp/9t,  is  defined  as 

^^p  _  P(.Z3)  Aj^  *  PB^S/i~AgP  I9B 

dt  mp 


where  m-  =  mass  of  the  projectile. 


Substituting  3Vp/9t  into  equation  (10) , 


P(z)  =  P 


Sr 


CAlvl  CAgAj^PUa) 


V^(zp)  V^(zp)m. 


/ 


Viz,  t) 
A(z,  t) 


dz 


CAjPr 


^b^ba^b  ^  _ 

V^(zp)mj,  V^(zp)m- 


rii^t^ldz 

J  Aiz.  t) 


CAbVpV^{z,  t) 
2V^  (zp)A^  (z,  t) 


(12) 
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Evaluating  equation  (12)  at  z  =  za  and  solving  for  P(za) 
and  defining  P(za")  =  P(za), 


.  CAlvi  v(z.  C)  g.  CA,A„P,  “ 

V^izp)  {  A(z,t)  ~  v^(za)in^- 

f  az 

A(z.t) 

CAgAj^  f  ^5^'  dz 
^*J  A(z,  t) 

1  +  ° 

( zp)  atp 

CAiPraa  7  Viz,  t)  VHza') 

V^(zp)mpJ  Aiz,t)  2V^  izp)  A^iza') 

1  +  ° 

V^izp)mp 

or 

P(za)  =  Z^O-PflB  *  ^ml  *  • 


where 


za’ 


f 


Viz,  t) 
hiz,  t) 


dz  , 


Oa  (za-) 


V^  jza') 
A^iza') 


^ao 


_ 1 _ 

^  ^  CAgAj^Q^iza-)' 
V^{zp)mp 
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= 


( zp) 


C>i(za‘)  + 


caIp, 


(zp)m. 


C?i(za')  - 


OAza-) 


2V^{zp) 


1  -*■ 


{ zp)  /n, 


Oi  (za') 


‘'a2 


(zp)nip 

^  ^  CA^j,0,(za-) 
( zp)  rrip 


Then 

P(z)  =1^3,  +[33(0  +  aJOPg  +  a5(C)Pa^]C>x(z)  +25(t){?2(^)  * 


where 


aj  ( t) 


CAlvl 

V^izp) 


^  CAb^zos 

V^{zp)/i7j,  v^izp)mp 


( zp)  nip  ( zp)  n7p 


ag  ( t) 


^g^A^aO 

V^{zp)nip 


bit) 


CAlvl  ^ 

2V^{zp) 


C?1  (z) 


r  y(z,  t) 
i  -fl(z,  t) 


3z  , 


C>2(z) 


V^iz) 

(z) 
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Now  the  area  from  the  aft  end  of  the  projectile  to  the  base  of  the 
projectile  will  be  examined. 


za  <  z  ^  zp 


Velocity  is  now 


u(z,  t) 


VpA^Viz.  t)  Aj,V^ 

V(zp)A(z,t)  A{z,t) 


(13) 


Taking  the  partial  derivative  of  u  with  respect  to  t  and  making 
some  substitutions 


dA(z,  t) 
dt 

dv(z,  t)  ^ 
dt 


=  0  . 


i 


dV(zp) 

dt 


= 


the  equation  for  3u(2,t)/9t  then  becomes 


duiz,  t)  _  VpA^Viz.t)  V^A^Aa 

dt  ~  V{zp)A  ^  V{zp)A{z,t) 
_  V^Alviz,  t)  _  Aj^Vp 
V^{zp)A{z,t)  A(z,t) 


Taking  the  partial  derivative  of  equation  ( 


j-j;  wxun  respect 


VA  fc) 

du(z,  t)  ^  dz _ 

dz 


V{zp)A{z,  t) 

_  ^pAgVjZft)  dA{z,t)  ^  Aj^Vp  dA[z,  t) 

V{zp)A^  dz  A^[z,t)  dz 

t) 
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du(z,  t)  ,  Vb  V{z,t)  aA(z,t) 

dz  ~  Vizp)  V(zp)  A^(z.t) 

+  3Alz,t) 

A^(z,t)  dz 


Since 


^  du(z,  t)  ,  ^..du(z,  t) 

<• — St~  — 5i~^ 


(14) 


equation  (14)  becomes 


+ 


c 

V^j,Viz,t) 

.  V^ApA^ 

Vizp) 

Vizp)Aiz.  t) 

Vizp)A_ 

C 

vlAlviz,  t) 

Aa 

V 

Vizp) 

izp)Aiz,  t) 

Aiz.  t) 

C  ( 

V^j,Viz,  t) 

]. 

vizp)  1 

Vizp)Aiz,  t) 

Aiz,  t)  } 

Vizp) 

c  ( 

VJipViz.  t) 

^kVp  ) 

V^j,Viz.  t)  dAiz,  t) 

Vizp)  1, 

Vizp)  Aiz,  t) 

Aiz,  t)  ) 

Vizp)A^ 

dz 

c 

(  V^gViz,t) 

^aVp 

.) 

^aVp 

dAiz,  t) 

Vizp) 

\  Vizp)  Aiz.  t) 

Aiz,  t) 

j 

A^iz,  t) 

dz 

I- 


(15) 


Solving  equation  (15)  for  8p  and  integrating, 


Z 


za* 


f  Viz,  t)  ^  CAJ^VJ,  r  dz 
V^(zp)  J  A{z,t)  Vizp)  J^A(z,t) 

^  CVpAe  r  V^jz,  t)  dA{z,  t) 

V^  izp)  J  ,A^(z,  t)  dz 

za 


2CVpAj^B  r  V(z,  t)  dA{z,t)  ^  CA^Vp  r  1  dA  (z,  t) 
V^izp)  J  A^iz.t)  dz  Vizp)  J  A^  iz,t)  dz 

(16) 
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Integrating  portions  of  equation  (16)  by  parts 


let 


V  =  I z,  t)  , 

du  =  , 

A^(z,  t) 


dv  =  2V(z,  t)A(z,  t)dz  , 


u 


1 

2A^(z,  t) 


then 


r  v^{z,  t) 
J  ,A^  (z,  t) 


dA(z,  t) 


V^(z,  t) 
2A^  (z,  t) 


Viz,  C) 
A(z,  t) 


dz 


V^(Z,  t)  ^  V^izan  ^  f  V(z,  t)  g, 
2A^iz,t)  2A^{za*)  J^^Aiz.t) 


and  letting 


V  =  V{z,  t)  , 

du  =  , 

A^iz,  t) 

dv  -  A(z,  t)dz  , 

u  - - ^ -  , 

2A^(z,  t) 


then 


r  v(z,  t) 

lA^iz.t) 


dA{z,  t) 


V(z, t)  1  +  f  dz 
2A^(z,t)za'  t) 
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v(z,  t)  ^  v(za*)  ^  r  dz 
2A^(z,t)  2A^(za^)  J^dA(z.C) 

JB^ 


Also, 


r  3a(z,  t)  _ _ 

,A^(z,  t)  2A^(z 


I -  \  =  - 1 _  .  _ 1 _ ] 

[z,  t),a‘  2\  A^(z,t)  A^{za*)j 


Therefore,  equation  (16)  becomes 


fap=  ..^Yp^b  f  yjz.  t)  3^  ^  f__dz 

J  V^(zp)  J  A(z,  t)  V(zd)  J  Alz. 


V(zp)  J  A(z,  t) 

za* 


CVpAa  V^{z,  t)  ^  V^izan 
V^(zp)  2A^(z,t)  2A^{za*) 


M  J  A(z,  t) 


2CV^Aj^g  Viz,  t)  ^  V(za*)  ^  If  clz,  J.. 

V^izp)  2A^(z,t)  2A^{za*)  2J  A 

L 

.  cAii^IiL  1  .  1  11 


V{zp)[2[  AUz,t)  A^iza*) 


Using  the  definition  for  9v(z,t)/3z,  and  P(za),  9V(z,t)/dz  becomes 

iV  _  ^A^aO^Br  *  ^A^al  *  ^A^az^B  ^B^BA  ~  ^B^ies 


fdP=- 


^^iAb  f  Viz,  t) 
v^izp)  J  ,Aiz,  t) 


dz  +  f 

Vizp)  J  Aiz,  t) 


CVlAl  vHz.t)  ^  CVlf^l  V^iza*) 
2V^izp)  AHz,t)  IV^izp)  A^iza*) 


+  f  V(g>  t)  0^  +  V(z,  £)  _  CVpAj^g  y(za*) 

izp)  j^,Aiz,  t)  V^izp)  A^iZft)  V^izp)  A^iza*) 

_  CVlA^g  r  dz  _  CA^Vl  ^  CAlvl 

1/2  /  J  Alz.  t)  OV(yr>'iA2  •?  1/f  ttiI  a 2  /  \ 


CVpA^B  f  dz  _  CAlvl  ^  CaIv^ 

V^izp)  J,A(z,C)  2V(zp)A^  *  2Vizp)AUza*) 
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Therefore,  P2(z)  is 


^2(2)  *  Piza*)  -  Q  (za^)  -  (^a*) 

vUzp)m^  '  vUzp)m^  ^ 


( zp)  m. 


0^(zan  - 


“p 

CA^baPb 

( zp)  m. 


OA  ^  P 

Q^iza^)  ^  QAza') 

v^{zp)m^^ 


-  V(2,t)  ^  CVpAp  v^iza*)  ^  CvlAl 

2V^  (zp)  A^(z,t)  (zp)  A*(za’)  v^  (zp)  ^ 

^  ^p^Ai^B  V(z,  t)  _  ^p^aAb  V(za*)  _  CVpA^iAj  r  3z 
V^(zp)  A^(z,t)  V^(zp)  A^(za*)  V^(zp)  J^^A(z,t) 

-  ^  cajvg 

2V(zp)A2(z,  t)  2V’(zp)A*(za*) 


Using  the  following  substitutions. 


Oi(za")  =  J 

zm* 


Viz,  C) 
A(z,  t) 


Bz  , 


Oj  (za*) 


/ 

za  * 


dz 


A(z,  t) 


the  equation  for  P2(z)  becomes 
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-f 


P^Kz)  =  P(za*) 

C^P Bz^3  ^  ^  '*’ 


+  a^lDOiiza*) 

+  CjiOojza*) 

CVpAaV^iza*) 

2V^(zp)A^(zan 


+ 


a«(t)P*C?i(^a*) 

+  C^{t)Pgpj{zan 

cvIaI 


2  (zp) 


Ooizp) 


+  l^(za*) 

V^(zp)  *  V^{zp)  A^{za*) 

2V(zp)A^iza*) 


CAa 

2V(zp) 


C>5  (zp) 


where 


03(0 


^A^al  _  ^jAa^zms 
v(zp)mp  V(zp)mp 


C^I^jAb 

( zp) 


^(t) 


^aAba  ^  ^A^a^ 

V(zp)mp  V{zp)mp 


C^it) 


CAa  Zap 

v{zp)mp 


Oiizp) 


1?4  ( zp) 


Os  ( ^P) 


V^(z,  t) 
A^{z.  t) 


Viz,  t) 
A^{z,C) 


A^iz.C) 


and  83,  a^,  a5,  are  the  same  as  stated  previously. 
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From  equations  (5)  and  (6) 


u(za*)  -  u(za') 


V^^Viza*) 
V{zp)  A(za*) 


A{za*) 


v^gViza-) 
V{zp)  A[za') 


letting 


V(za*)  =  V{za')  =  v(za) 


^  V^^V(za)f  1 _ _  A^v^ 

V{zp)  A(za-)j  A(za*)  ' 


and  noting 


= 


A{za')  -  A(za*) 
A  ( za  * )  A  ( za  * ) 


With  the  Kooker  anaylsis, 

fk  =  jump  in  pressure  across  the  boattail. 


An  analysis  by  Kooker  (April  1991)  indicates  that  the  pressure 
drop  across  the  boattail  (which  in  a  one-dimensional  analysis  is 
equivalent  to  determining  the  pressure  drop  across  a  moving 
discontinuity  in  area,  Figure  A-2)  is  given  by: 
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Mass  balance; 


Px(“l  -  Vp)Al  =  P2(U2  -  Vp)A2  •  A 


(17) 


Momentum  balance: 


P2^2(*^2  "  ~  Pl^l(*^l  “  ~  ~  ~ 


=  -(Pj  -  p,) 


Aj  +  Aj 


(18) 


since 


=  i{Pl  +  i’2)- 


Thus, 


{^2  -  Pi) 


2p(Ui  -  Vg)A^{U2  -  Ui) 

Ai  +  Aj 


(19) 


where 

Ai  =  A(za-) , 

Aj  =  A(za*) , 

=  u(za') , 

Uj  =  u{za*) , 

or 

P^Ua*)  -  Pi(za-)  =  f ic  = 


-2CVlA{za-)  Aj^ 

AgV(za)  Y 

V(zp)  A(za*) 

V{zp)A(za~) 

A{za')  +  A{za*) 


(20) 


Therefore,  P2(z)  is 
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p^iz)  =  P^(za*)  +  P^(za-)  *  fk 


P^iz)  =  z.oPa,  +  +  f Ic  +  a^O^{.za*) 

+  a4Pai?i(za*)  +  ajPg^OiCza*)  +  C^O^iza*) 

*  C^PaO^iza*)  +  C^Pa^Ojiza*)  +  ^(OOjCza*)  +  AiCi4(za*) 

+  Ji05(za*)  +  ;ci  . 


Let 

Pg  =  PjCZp) 


^h(1  -  ^*2  -  a4C)i(zp)  -  C^Oj(zp))  = 

■Par(^«o  +  a5(?i{zp)  +  C^Oiizp))  + 

+  fk  +  a^^Oiizp)  +  CjOjCzp)  +  l)(  t)02(  zp)  +  *i04(zp) 

+  JiOsizp)  *  K 


•Pfl-^2  ~  ^3^BZ  *  • 


where 

ll  =  Zai  *  *  ^iOx(zp)  +  C3C>3(zp) 

+  Jb(t)(?2(zp)  +  h^o^izp)  +  JiOs^zp)  *  ^1  . 


•^2  ~  1  ~  ^a2  ^40^  izp)  ~  ^4l?3  i  Zp)  , 

I3  =  z^o  *  a5C>i(zp)  +  C5p3(zp)  , 


Olizp) 


f 

WM* 


V{z,  t) 
A(z,  t) 


Bz  , 
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-CAj,VB^V(za)Aj^  CV^Aiv^iza) 

V^(zp)A^{za)  2V(zp)A^(za)  2V^  (zp) A^  (za) 

_  CAj^^p 

vMzp)  ' 

.  _  -CA^V^ 

2V(zp) 


We  are  now  at  a  point  where  the  projectile  base  and  breech 
pressure  can  be  determined.  For  use  in  lumped-parameter  interior 
ballistic  models,  the  gradient  equation  is  usually  cast  in  terms  of 
the  mean  pressure  (P^j)  : 
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Substituting  the  pressure  distribution, 


= 


0 


jA(z,c)Pgj.*  J  A(z.  c)  a^ic)  OiBz  j  A(z,  t)  a^(c)  O^Pg^dz 


v(zp) 


0  Q  zp 

j  A{z,  t)  a^{t)  Oy^P^z  *  ^  A{z,  t)b(t)0^dz  *  j  A(z ,  t)  z^^Pg^dz 


za' 


za^ 


V(zp) 

IP  IP 

fA(z.  t)z^,dz  +  I  /  A(z,  t)  fkdz 


IP 


za* 


za 


Vizp) 


Zp  zp 

j  A(z,  t)  a^Oi  (za'^)dz  +  J A(z,  t)  a.^P^y  {za*)dz 


za* 


za* 


V(zp) 


IP  Ip 

j  A(z,  t)  a^PgrOi  (za*)dz  +  J  A(z,  t)  C3O3  (za*)dz 


za* 


za* 


V(zp) 


tp  Ip 

f  A(z,  C)C^P^J(zan^z  *  J  AKz,  t)  C^Pg^O^iza^  dz 


V(zp) 


Ip  zp 

J  A(z,  t)  b{t)  O^dz  +  J  A(z,  t)  hyO^dz 


za* 


za* 


V{zp) 

zp  zp 

j  A(z,  t)  jyO^dz  +  jAiz.t)k^dz 


za* 


za* 


v(zp) 


Rearranging  and  making  some  substitutions,  the  equation  becomes 
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^  _  33(007  >  t{t)0^  +  z^i(v(zp)  -  v(za)) 

“  Vizp) 

a^{t)0^  +  fkiv(zp)  -  v(za))  +  C^iOQg 
*  vizp) 

h^O^iza*)  +  j^O^iza)  +  k^jvjzp)  -  v(za)) 

Vizp) 


+  P, 


viza)  +  35(007  +  a5(t)09 


+  P« 


Vizp) 

Zao(V(zp)  -  V(z3)  )  +  C^it)Qg 


+  p, 


vizp) 

a^it)Oj  *  a^it)Q^  +  z^2  ^  ^P'>  ~  v(za))  C^it)Q^ 

vizp) 


where 


v^iz,  O 
Aiz,  t) 


dz  , 


0  0 


Oa  =  ^Aiz,  t)  f  . 

za'  za  * 


O9  -  "[Aiz,  t)  f  dzdz  , 

za'  za' 


and  Q2/  Q4,  and  Q5  are  as  stated  previously. 
Or 
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where 


a^[t)Q,{za)  *a^{t)Q^[zp)  -»•  i3(  t)  (zp) 

vTzp) 

^  z^^ivjzp)  -  V(za))  *  fk{v(zp)  -  v{za)) 

Vizp) 

+  ^  ^lOi  ( jiOj  ( ^P) 

Vizp) 

+  ^1  ~  ) 

V(zp) 


_  ^4  ( 0?  (ga)  •*•  ( t)  gg  (zp) 

^  v(zp) 

z^3(v(zp)  -  Vjza))  ■*■  C^it)  Qf^izp) 
Vizp) 


_  -*•  ag  ( t)  ^7  (za)  *  a^Q^izp) 

3  vCzp) 

j  z^o^Vizp)  -  Viza))  •>•  CjOet^^p) 
\^{zp) 

Therefore, 


Subtracting  the  two  previous  equations  and  solving  for  Pg, 
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Pb- 


The  energy  of  the  fluid  is  represented  by 

dE  -  ~u^  (z,  t)  dm  , 
dm  =  pA{z,  t)dz  . 


Integrating  from  0  to  zp, 


(z,  t) 
2 


pA(z,  t)  dz 


za"  zp 

=  J  (z,  t)  A  (z,  t)dz  *  j  (z,  C)A(z,  C)3z 

0  2a* 


C 

2V{zp) 


J  AgVpV^(Z,t) 
I  A(z,  t)  (zp) 


+  C  J(  VpAsV(z,t) 
2V{zp)  J  [v{zp)A(.z.  t) 


A[z,  t) 


j  A(z,  t)0z 


f  V^iz,  t)  ^ 
2V^{zp)  {  A(z,  t) 


CVp  J  A^v^ 
2V{zp)  izp)A(z,  t) 


+ 


cvl  J{  2Aj^AbV{z.  t) 
2V(zp)  i  i,  V{zp}A[z,t) 

9M  * 


+ 


B  2 


A(z,  t) 


dz 


49 


_  r  v^iz,  t)  ^  CVpAe  ?  (z,  t)  g, 

2V^(zp)  J  A(z.t)  2V^(zp)  J  Aiz.t) 

u  za 

_  CVpAaAb  7  t)  g_  ^  7 

V^(zp)  J  Aiz.t)  2Vizp)  J  Aiz.t) 

za  *  za  * 


Substituting  Q2(zp),  Q3(zp),  and  Qg  (which  were  previously  defined) 
into  the  energy  of  fluid  equation, 


2  ( zp) 


Ce 


i  zp) 


0i  i  zp) 


+ 


CV^aI 
2Vi  zp) 


O^izp)  . 


Volume  and  area  terms  for  a  straight  tube  (Lagrange  tube) : 

Viza)  =  Agza 

ViZP)  =  AgZa  *  (Zp-Za)  Ajgyj 


za‘ 


P  _  7  t)  _  (^B^a  +  (^P  -  ^a)A3*)2  (A^za) 


13  (za)  =  =  za^ 

^  A2(za) 


C?2(2p)  = 


V^(zp)  ^  (^jza  +  Ag^izp  -  za))^ 
A^izp) 


A  ^ 
■"BA 


v(zp) 

z  2 
•AflA 
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O^(zp) 


Oi(zp)  = 


‘  dz 

_  zp  -  za 

A(z,  t) 

1 

Agj^ 

V(zp) 

_  V{zp) 

A^izp) 

a  2 

■Aba 

1 

1 

A^(z,t)  aI 


Agza^  (Agza  +  (za  -  za)A^^ 
-  + 


3Ai, 


3A|» 


zd 


ft(za)  =  /A(z,t)/||f^a^az. 


Agza^ 


za  *  za  * 


0,(zp)  =  f  A(z.  t,  f  {^ez3z 
za^  za' 

{AgZa  +  Aj^(zp  -  za)f  _ 


-  -  i52ffl!(zp-  za) 


®-"aA 


2A 


BA 


The  energy  term  for  the  Lagrange  tube; 


za)Agj^'^ 
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dE  =  (z,  t)  dm 

A 


dm  =  pA(z,  t)dz 


gp  za  zp 

jdE=  ^  Ju^  (z,  t)  A(z,  t)dz  +  ^ju^(z.t)A(z,t)dz 

0  0  za 


-  -gf  Aiz.t)a. 

2  J  Vizp)  ^A^  (z,  t) 

J(  _  _A^Y 

J\v(zp)A(z,t)  A(z.t))  ^  ’ 


za 


Cfljvg  Jz^A^ 

2V^(zp)l  Ai 

(z  -  za)AgJ^) 


t)dz 


dz 


2V(zp) 


J(  f 

za' 


Vizp) 


CAjv^  za3 


-  Aj,V^ 


2V^{zp)  3 


v*(zp)  ■  v[zp)  * 
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2  ( zp)  3 


C 

2V(zp)Agf^  V^(zp) 


j  (Agza  +  (z  -  za)AgJ^)^^z 


C  2A.Vj?A«?  . 

2V'(zp)A„  VUp)  /*-*•”  * 


iv^/az 


Substituting  Qi(zp),  Q3(zp),  and  QeCzp)  into  the  equation,  the 
energy  equation  becomes 

JdE  =  -  ^'0,(zp)  .  . 

<  2  { zp)  ( zrt)  2  V ( zp) 
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Intentionally  left  blank. 
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APPENDIX  B: 

USER’S  MANUAL  AND  CODE  LISTING  FOR  RGA 
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Intentionally  left  blank. 
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o  o 


program  ibrga 

common  nsll,  kpr,  fracsl(lO),  dsdxsKlO),  surfsl(lO),  nslp(lO), 

1  tsl(lO),  plsrch,  phase,  pmean,  bbr(lO),  abr(lO),  deltat,  y(20), 

1  igrad 

character  cutfil  *  10,  bdfile  *  10,  style  *  10 
character  title  (15)  *  4,  vsn  *  4 

dimension  br{10),  trav{10),  rp(lO),  tr(lO),  forcpdO),  tempp(lO), 
1  covp(  10) 

dimension  chwp(lO),  rhopilO),  gamap{10),  nperfs(lO),  glenp(lO), 

1  pdpdO  ),  gdiapdO),  alphadO,  10),  beta  (10,  10),  pres  (10,  10) 
dimension  a(4),  b(4),  a]c(4),  d(20),  p(20),  z(20),  frac(lO), 

1  surf(10  ),  volp(lO),  dsdx(lO),  nbr(lO),  ibo(lO),  tbo(lO), 

1  d2xdt2  (10)  ,  tngdO) 

real  lambda,  jlzp,  j2zp,  j3zp,  j4zp,  jlzb,  j2zb,  j3zb,  j4zb 
real  11,12,13 

dimension  chdist(6),  chdiam(6),  bint (10),  projtr(20),  projms(20) 
dimension  nsl(lO),  surfodO),  dsdxn(lO) 
data  pi/3 . 14159/vsn/' 4hboat ' / 
c 
c 

c  USER'S  MANUAL  FOR  IBRGA 

c 

c 

c  IBRGA  relies  on  an  input  database  consisting  of  all 

c  numerical  parameters  essential  for  running  the  code.  Values  may 
c  be  in  metric  units  or  in  Imperial  units,  but  must  be  consistent 

c  throughout  a  dataset.  Below  is  a  compilation  of  a  typical  data 

c  base  showing  the  name  and  location  of  each  parameter.  The  names 

c  for  the  numerical  values  are  prefixed  with  an  alphabetical 

c  designator  corresponding  to  the  position  at  which  the  data  is  to 
c  appear,  that  is,  from  left  to  right.  The  data  may  be  separated 
c  by  blanks  or  commas.  Measurement  units,  if  any,  are  shown  to  the 

c  right  of  each  input.  In  general,  metric  units  of  weight  and  mass 

c  are  the  meter  and  kilogram,  respectively;  corresponding  Imperial 
c  units  are  the  inch  and  pound.  The  only  exceptions  are  Imperial 

c  units  of  propellant  impetus,  which  are  foot-pounds  per  pound  mass, 

c 
c 

c  title  card  -  up  to  60  characters  of  title  and  identification 

parameter  information  and  placement: 
c 

cABCDEFGHIJK 

c 

c  record  1  Metric  Imperial 

c  A.  -  chamber  volume  (when  canister  (m**3)  (in**3) 

c  model  is  used,  this  will  be  the 

c  volume  after  canister  bursts) 

c  B.  -  groove  diameter  (m)  (in) 

c  C.  -  land  diameter  (m)  (in) 

c  D.  -  groove/land  ratio  (land,  groove, 

c  and  groove/land  ratio  used  to 

c  calculate  the  tube  bore  area) 

c  E.  -  twist  (units  are  turns/caliber) 

c  F.  -  projectile  travel  (m)  (in) 

c  G.  -  gradient  switch  (integer  value 

c  designating  the  gradient  equation 

c  (1  =  Lagrange,  2  =  Chambrage, 
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o  o 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


3  =  Two-phase,  4  =  RGA, 

5  =  Lagrange  w/  bt,  6  =  Cham,  w/  bt) 

H.  -  variable  projectile  mass  switch 

(0=no,  l=yes) 

I .  -  igniter  canister  model  switch 

(0=no,  l=yes) 

J.  -  friction  factor  (normally  1  for 

granular,  0.01  for  stick  and 
0.1  for  partially  cut  propellant; 
only  used  when  gradient  •  3  or  4) 

record  la  (read  if  and  only  if  gradient  »  5  or  6) 

A.  -  boattail  diameter  (m)  (in) 

B.  -  boattail  length  (m)  (in) 


record 

A.  - 

B.  - 

C.  - 


lb  (Read  if  and  only  if  gradient  •  2  or  4  or  6) 


number  of  point  pairs  to  describe 
chamber  geometry,  integer  I  <=>  5 
initial  distance  from  breech 

(m) 

(in) 

(must  be  0.0) 

diameter  at  initial  distance 

(m) 

(in) 

c 

c 

c 


c 

, 

- 

Ith  distance  from  breech 

(m) 

(in) 

c 

(initial  position  of  the  base 

c 

of  the  projectile) 

c 

- 

Ith  diameter  at  Ith  distance 

(m) 

(in) 

c 

(used  to  calculate  bore  area  - 

c 

overrides  record  1  groove  and 

c 

land  diameter  specifications) 

c 

(Note:  chamber  geometry  is  used 

c 

to  calculate  the  chamber  volume 

c 

which  overrides  record  1  chamber 

c 

c 

volume  description.) 

c 

c 

record 

2 

c 

A. 

- 

projectile  mass 

(kg) 

(lb) 

c 

B. 

- 

switch  to  calculate  energy  lost 

c 

to  air  resistance,  an  integer 

c 

either  0  =  no  loss,  or  1  =  loss 

c 

C. 

- 

fraction  of  bore  resistance  work 

c 

used  to  heat  tube  (0 . 0<=f<=l . 0) 

c 

D. 

— 

gas  pressure  ahead  of  projectile 

(MPa) 

(psi 

c 

c 

record 

2A  (Read  if  and  only  if  variable 

projectile 

mass 

c 

switch  is  1) 

c 

A. 

- 

number  of  point  pairs  to  describe 

c 

variable  projectile  mass  =<  20 

c 

B. 

- 

initial  projectile  travel 

(m) 

(in) 

c 

(conceptually  should  be  0.0) 

c 

C. 

- 

initial  projectile  mass 

(kg) 

(lb) 

c 

(overrides  value  from  record  2) 

c 

D. 

- 

projectile  travel  at  which  first 

(m) 

(in) 

c 

mass  change  occurs 

c 

E. 

- 

new  projectile  mass  at  travel  D. 

(kg) 

(lb) 

c 
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c 


c 

X  - 

i-th  projectile  travel  where  mass 

(m) 

(in) 

c 

change  occurs 

c 

y  - 

i-th  new  projectile  mass  value 

(kg) 

(lb) 

c 

c 

record 

2B  (Read  if  and  only  if  igniter 

c 

canister  model  switcn  =  1) 

c 

A.  - 

pressure  at  which  the  igniter 

(MPa) 

(psi) 

c 

canister  will  burst 

c 

B.  - 

volume  of  igniter  canister 

(m**3) 

(in**3) 

c 

(used  as  chamber  volume  until 

c 

burst  pressure  achieved) 

(in) 

c 

C.  - 

canister  diameter  (assumes  a 

(m) 

c 

right  circular  cylinder) 

c 

c 

record 

3 

c 

A.  - 

number  of  pairs  of  barrel 

c 

resistance  points  (integer  <=  10) 

c 

B.  - 

bore  resistance 

(MPa) 

(psi) 

c 

c 

C.  - 

travel 

(m) 

(in) 

c 

c 


c 

- 

Jth  bore  resistance 

(MPa) 

(psi) 

c 

• 

— 

Jth  travel 

(m) 

(in) 

c 

c 

record 

4 

c 

A. 

- 

mass  of  recoiling  parts 

(kg) 

(lb) 

c 

B. 

- 

number  of  recoil  point  pairs 

c 

(must  be  an  integer  =  2) 

c 

C. 

- 

recoil  force  (force  to  overcome 

(N) 

(lb) 

c 

before  recoil  start  -  rod  preload) 

c 

D, 

- 

time  of  rod  preload  (must  be  0.0) 

(s) 

(s) 

c 

E. 

- 

recoil  force  (constant  resistive 

(N) 

(lb) 

c 

force  after  rise  time) 

c 

F  , 

- 

rise  time  (time  to  go  from  recoil 

(s) 

(s) 

c 

start  to  constant  resistive 

c 

recoil  force) 

c 

c 

record 

5 

c 

A. 

- 

free  convective  heat  transfer 

(W/m**2/K) 

( in-lb/ in*  *2 

c 

coefficient 

/s/K) 

c 

B. 

- 

chamber  wall  thickness  (wall 

(m) 

(in) 

c 

depth  whic).  is  heated  uniformly) 

c 

C. 

- 

heat  capacity  of  chamber  wall 

(J/kg/K) 

(in-lb/lb/K) 

c 

D. 

- 

initial  temperature  of  tube  and 

(K) 

(K) 

c 

chamber  walls 

c 

E . 

- 

heat  loss  coefficient  (usually  1., 

c 

but  may  be  set  to  0.0  in  order  to 

c 

eliminate  heat  Joss) 

c 

F. 

- 

density  of  chamber  wal’ 

(kg/m**3) 

(lb/in**3) 

c 

c 

record 

6 

c 

A. 

- 

impetus  of  igniter 

(J/kg) 

(ft-lb/lb) 

c 

B. 

- 

adiabatic  flame  temperature  of 

(K) 

(K) 

c 

igniter  material 

c 

C. 

- 

covolume  of  igniter 

(m**3/kg) 

{ in*  *3/lb) 

c 

D  . 

- 

ratio  of  specific  heats  of  igniter 
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(lb) 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c* 


E.  -  mass  of  igniter  {)cg) 


record 

7 

A.  - 

number  of  propellants 

(integer  <=  10) 

record 

8 

A.  - 

impetus  of  propellant 

(J/kg) 

(ft-lb/lb) 

B.  - 

adiabatic  flame  temperature 

(K) 

(K) 

C.  - 

covolume  of  propellant 

(m**3/kg) 

(in**3/lb) 

D.  - 

ratio  of  specific  heats 

E.  - 

mass  of  propellant 

kg) 

(lb) 

F.  - 

density  of  propellant 

(kg/m*  *3 ) 

(lb/in**3) 

G.  - 

propellant  form  function  indicator 

(integer;  may  be  one  of: 

0  solid  cylindrical  grain 

1  single-perf  cylindrical  grain 

2  spherical  grain 

7  seven-perf  cylindrical  grain 

15  nineteen-perf  hexagonal  grain 

19  nineteen-perf  cylindrical  grain) 

H  .  - 

length  of  propellant  grain 

(m) 

(in) 

I.  - 

diameter  of  perforations  in  the 

(m) 

(in) 

propellant  grains  (ignored  if  not 
required,  but  must  be  present) 

J.  - 

outside  diameter  of  propellant 

(m) 

(in) 

grain  (for  the  hexagonal  grain 
it  is  the  distance  between 
rounded  corners) 

(Record  8  repeated  for  each  propellant) 


record  9 


A. 

-  number  of  burning  rate  triplet 

points 

(integer  J  <=  10) 

B. 

-  exponent 

C. 

-  coefficient 

(m/s-MPa**e) 

(in/ s-psi**e) 

D. 

-  pressure  (upper  pressure  limit 

(MPa) 

(psi) 

for  which  the  previous  exponent 
and  coefficient  are  valid) 


-  Jth  exponent 

-  Jth  coefficient 

-  Jth  pressure  (if  pressure  should 

exceed  this  limit,  then  this 
burning  rate  equation  is  used 
for  all  higher  pressures) 


(m/s-MPa**e)  (in/s-psi**e) 
(MPa)  (psi) 


(Record  9  repeated  for  each  propellant) 


record  10 

A.  -  integration  time  increment  (ms) 

B.  -  print  increment  (ms) 

C.  -  upper  limit  on  integration  time  (ms) 

for  stopping  calculation 


(ms) 

(ms) 

(ms) 
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c*  conversion  factors  for  imperial  units  =>  metric  units 
c* 


c 

length 

inches 

*  .0254 

=> 

meters 

c 

mass 

lb  *  . 

45359237 

=> 

kilograms 

c 

area 

in''2  *  . 

00064516 

=> 

m^2 

c 

volume 

in"3  *  0 

1 .000016387064 

=> 

m'^S 

c 

pressure 

lb/in"2 

*  6894.757 

=> 

pascals 

c 

velocity 

ft/s  / 

3.28083 

=> 

m/s 

c 

in/s  * 

.0254 

=> 

m/s 

c 

energy 

ft-lb  * 

1.3558179 

=> 

joules 

c 

in-lb  * 

0.1129848 

=> 

joules 

c 

density 

lb/in''3 

*  27679.9 

=> 

kg/m''3 

c 

force/mass 

(ft-lb)/lb  *  2.989067 

=> 

j/kg 

c 

CO volume 

in^3/lb 

/  27679.9 

=> 

m''3/kg 

c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 

c 

c 


source:  engineering  design  handbook  metric  conversion  guide 
darcom  pamphlet  706-470,  july  1976 

density  grams/cc  *  1000.  =>  kg/m'^S 

call  gettim (ihr , imin, isec, ihuns) 

write(  *,  830) 
read(  *,  840)bdfile 

open(unit  =  2,  err  =  810,  file  =  bdfile,  status  =  'old',  iostat 
1  ios) 
nzp=0 
rewind  2 
write  (  *,  850) 
read(  *,  840)outfil 

open(unit  =  6,  err  =  820,  file  =  outfil) 
dol0i=l,  20 


p(20) 

=  0. 

y(20) 

=  0. 

z  (20) 

=  0. 

d(20) 

=  0. 

10  continue 

write  (  *, 

870) 

read(  *, 

840  ) 

mode  =  0 

if (style  (1 : 1)  .eg. 'm'  .or.  style (1 : 1) .eq. 'M' )  mode  =  1 

if (style (1 : 1) .eg. 'e'  .or.  style (1 : 1 ). eq. ' E' )  mode  =  2 

if (mode . eq . 0)  write (  *,  880) 

if (mode .eg. 0)  stop 

read(2,  885)title 

write(6,  1236) title, vsn 

write(6,  860)bdfile 

read(2,  end  =  790,  err  =  800)cham,  grve,  aland,  glr,  twst, 

1  travp,  igrad,  ivpm,  ihl,  fsO 
if (igrad.gt . 1) go  to  20 
write(6,  890) 
igrad  =  1 
go  to  140 

define  chambrage  assumes  nchpts=number  of  points  to  define 
chamber  >  or  =  2  <  or  =  5  (?) , chdiam(i)  defines  chamber  diameter 
at  chdist  (i)  chamber  distance.  chdiam (nchpts)  is  assumed  to  be 
the  bore  diameter  and  chdist (i)  is  assumed  to  be  0,  i.e.  at  the 
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c 

c 


breech,  assumes  truncated  cones. 


c 

c 

c 


c 

c 


20  if ( igrad. eq . 3) go  to  130 
i f ( igrad . eq . 4 ) go  to  30 
i f ( igrad . GE . 5) go  to  25 
write(6,  900,  err  •=  800) 
go  to  40 

25  read  (2  ,  *,  end  =  790, err  =  800)  btdia,btlen 
i f (mode . eq . 1 ) then 
btvol=pi*btdia’*btdia/4  .  *btlen 
write (6, 955) btdia, btlen, btvol 

955  f ormat ( / , lx, ' bDattail  diameter  m  ' ,  el6 . 6/lx,  ' boattail  length 
&m  ', el6 . 6/lx, ' boattail  volume  m**3  ',el6.6,/) 
else 

btvol=pi*btdia*btdia/4 . *btlen 
write(6, 965, err=800) btdia, btlen, btvol 

965  f ormat {/, lx, ' boat  tail  diameter  in  el6 . 6/ lx,  '  boattail  length 
sin  ', el6 . 6/lx, ' boattail  volume  in**3  ',el6.6,/) 
btdia=btdia*0 .0254 
btlen=btlen*0 .0254 
btvol =btvol*0 .0254**3 
cham=cham*l . 6387064e-5 
endif 

35  if (igrad. eq. 5) then 
write (6, 975) 

975  format ( lx, ' using  lagrange  with  boattail  gradient') 
go  to  140 
endif 

if (igrad. eq . 6) then 
write (6, 995) 

995  format ( lx, ' using  chambrage  with  boattail  gradient') 
go  to  40 
endif 

30  write(6,  910) 
go  to  40 

40  read(2,  *,  end  -  790,  err  -  800)nchpts,  (chdist(i),  chdiam(i),  i 
1=1,  nchpts) 
if (mode .eq. 1) then 

write(6,  920,  err  =  800) (chdist (i) ,  chdiam(i),  i  =  1,  nchpts) 
goto  60 

else 

write(6,  925,  err  -  800) (chdist (i) ,  chdiam(i),  i  -  1,  nchpts) 
do  50  i  =  1,  nchpts 

chdist(i)  =  chdist(i)  *  0.0254 
chdiam(i)  =  chdiam(i)  *  0.0254 

50  continue 

endif 

calculate  chamber  integrals  and  volume 

60  if (nchpts . gt . 5)  write(6,  930,  err  -  800) 
if (nchpts .gt . 5) nchpts  -  5 
bore  =  chdiam (nchpts) 

if (chdi St  ( 1 ) . ne . 0 . 0 ) write  (6,  940,  err  =  800) 
chdist (1)  =0.0 
do  54  1=1, nchpts 
chdist(I)=0. 01*chdist (I ) 
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c54  chdiam ( I ) =0 . 01*chdiam ( I ) 

c  calculate  chamber  integrals  and  volume 

if (nchpts . gt . 5)  write (6, 44, err=30) 

44  f ormat ( lx, ' use  first  5  points') 

if (nchpts .gt . 5) nchpts=5 
bore=chdiam (nchpts) 

if (chdist (1) . ne . 0 . 0) write (6, 45, err=30) 

45  format (lx,'  #  points  ?  ') 

chdist (1) =0 .0 
ptl=chdist  (nchpts) 
btd=btdia 

btl=btlen 

call  j int (btd , bt 1 , pt 1 , pt 1 , nchpt s , chdi st , chdiam, bint , bvol ) 

41  cham=bvol+btvol 

write (6,47,err=30)bint(l) ,bint (3) ,bint (4) 

format (lx, 'bint  1  =  ',el4.6,'  bint  3  =  ',el4.6,'  bint  4  =  ',el4. 
&6) 

chmlen=chdist (nchpts) 
go  to  140 
130  write  (6,  950) 

140  if (mode .eq. 1) then 

write  (6,  960,  err  =*  800)cham,  grve,  aland,  glr, 

&  twst,  travp,  igrad,  ivpm,  ihl,  fsO 
cham=cham-btvol 
endif 

if (mode .eq. 2) then 
cham=cham/l . 6387064e-5 

write  (6,  970,  err  =  800)cham,  grve,  aland,  glr,  twst,  travp, 

1  igrad,  ivpm,  ihl,  fsO 

cham  =  cham  *  1.6387064e  -  5 
cham=cham-btvol 
grve  =  grve  *  0.0254 
aland  =  aland  *  0.0254 
travp  =  travp  *  0.0254 
endif 

read(2,  *,  end  =  790,  err  =  800)prwt0,  iair,  htfr,  pgasO 
if (mode.eq. 1) then 
prwt  =  prwtO 
pgas  =  pgasO  *  1.0e6 
el seif (mode . eq. 2) then 

prwt  =  prwtO  *  0.45359237 
pgas  =  pgasO  *  6894.757 
endif 

if (ivpm.eq, 1) then 

read(2,  *  ) nvpmp,  (projtr(i),  projms(i),  i  =  1,  nvpmp) 
if (mode .eq. 1) write (6,  980)nvpmp,  (projtr(i), 

1  projms(i),  i  =  1,  nvpmp) 

if (mode . eq . 2 ) then 

write  (6,  985) nvpmp,  (projtr(i),  projms(i),  i  =  1,  nvpmp) 
do  150  i  =  1,  nvpmp 

projtr(i)  =  projtr(i)  *  0.0254 
projms(i)  =  prcjms(i)  *  0.45359237 
150  continue 

endif 

prwt  =  projms(l) 
prwtO  =  prwt 

if (mode .eq. 2)  prwtO  =  prwt  /  0.45359237 
endif 


write(6,  1050) 
if  { ihl . eq . 1 ) then 

read (2,  *  )burstp,  highv,  highd 

if (mode .eg. 1) write (6,  990)burstp,  highv,  highd 

if (mode .eq. 2) then 

write  (6,  1000)burstp,  highv,  highd 
burstp  =  burstp  *  0.006894757 
highv  =  highv  *  1.6387064e  -  5 
highd  =  highd  *  0.0254 
endif 

burstp  =  burstp  *  l.e6 

areaw  =  4 .  *  highv  /  highd  +  pi  *  highd  *  highd  /  2 . 
endif 

read(2,  *,  end  =  790,  err  =  800)npts,  (br(i),  trav(i),  i  =  l,npts) 

read(2,  *,  end  =  790,  err  =  800)rcwt,  nrp,  (rp(i),  tr(i),  i=l,nrp) 

read(2,  *,  end  =  790,  err  =  800)ho,  tshl,  cshl,  twal,  hi,  rhocs 

read(2,  *,  end  =  790,  err  =  800)forcig,  tempi,  covi,  gamai,  chwi 

read (2,  *  ) nprop 

read(2,  *,  end  =  790,  err  =  800) (forcp (i) ,  tempp(i),  covp(i), 

1  gamap(i),  chwp(i),  rhop(i),  nperfs(i),  glenp(i),  pdp(i),  gdiap(i) 
1  ,  i=  1,  nprop) 
if (mode . eq . 1 ) then 

write  (6,  1010,  err  =  800)prwt0,  iair,  htfr,  pgasO 
write(6,  1030,  err  =  800)npts,  {br(i),  trav(i),  i  =  1,  npts) 
write(6,  1060,  err  =  800)rcwt,  nrp,  (rp(i),  tr(i),  i  =  1,  nrp) 
write  (6,  1080,  err  =  800)ho,  tshl,  cshl,  twal,  hi,  rhocs 
write  (6,  1100,  err  -  800)forcig,  tempi,  covi,  gaunai,  chwi 
write(6,  1236) title, vsn 
write  (6,  1120) nprop 

write(6,  1130,  err  =  800)  (i,  forcp(i),  tempp(i),  covp(i), 

1  gamap(i),  chwp(i),  rhop(i),  nperfs(i),  glenp(i),  pdp(i), 

1  gdiap(i),  i  =  1,  nprop) 
endif 

if (mode .eq. 2) then 

write  (6,  1020,  err  =  800)prwt0,  iair,  htfr,  pgasO 
write(6,  1040,  err  =  800)npts,  (br(i),  trav(i),  i  =  1,  npts) 
write(6,  1070,  err  =  800)rcwt,  nrp,  (rp(i),  tr(i),  i  =  1,  nrp) 
write  (6,  1090,  err  =  800)ho,  tshl,  cshl,  twal,  hi,  rhocs 
write  (6,  1110,  err  =  800)forcig,  tempi,  covi,  gamai,  chwi 
write(6,  1236) title, vsn 
write (6,  1120) nprop 

write  (6,  1140,  err  =  800)  (i,  forcp(i),  tempp(i),  covp(i), 

1  gamap(i),  chwp(i),  rhop(i),  nperfs(i),  glenp(i),  pdp(i), 

1  gdiap(i),  i  =  1,  nprop) 
endif 

do  170  j  =  1,  nprop 

read(2,  *,end=  790,  err  =  800)nbr(j),  (alpha(j,  i) ,  beta(j,i) 
1  ,  pres(j,  i),  i  =  1,  nbr(j)) 

if (mode .eq. 1) write (6,  1160)nbr( j) 
if (mode .eq. 2) write (6,  1170) nbr ( j) 
do  160  i  =  1,  nbr(j) 

if (mode .eq. 1) write (6,  1180)  alpha(j,  i) ,  beta(j,  i) , 

1  pres(j,  i) 

if (mode .eq. 2) then 

write(6,  1180)  alpha(j,  i),  beta(j,  i),  pres(j,  i) 
rate  =  beta(j,  i)  *  pres(j,  i)  **  alpha (j,  i) 
pres(j,  i)  =  pres(j,  i)  *  0.006894757 
beta(j,  i)  =  0.0254  *  rate  /  pres(j,  i)  **  alpha(j,i) 
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endif 

160  continue 

170  continue 


convert  units  to  program  requirements 

do  180  i  =  1,  npts 

if (mode.eq. l)br (i)  =  br(i)  *  l.e6 
if (mode .eq. 2) then 

br(i)  =  br(i)  *  6894.757 
trav(i)  “  trav(i)  *  0.0254 
endif 

180  continue 

do  200  j  =  1,  nprop 
if (mode . eq . 2 ) then 

forcp(j)  =  forcp(j)  *  2.989067 
covp(j)  =  covp(j)  /  27679.9 
chwp(j)  =  chwp(j)  *  0.45359237 
rhop(j)  =  rhop(j)  *  27679.9 
glenp(j)  =  glenp(j)  *  0.0254 
pdp(j)  ■=  pdp(j)  *  0.0254 
gdiap(j)  =  gdiap(j)  *  0.0254 
endif 

do  190  i  =  1,  nbr(j) 

pres(j,  i)  =  pres(j,  i)  *  l.e6 
190  continue 

200  continue 

if (mode . eq. 2) then 

do  210  i  =  1,  nrp 

rp(i)  =  rp(i)  0.1129848 
tr(i)  =  tr(i)  *  0.0254 
210  continue 

c  conversion  factor  for  free  convective  heat  transfer  coeff 
c  w/m**2  *  (0.00064516  m**2/in**2)  *  (1.0/1.3558179  ft-lb-s/w)  * 
c  (12.0  in/ft)  =  0.005710147  in-lb/in**2-s 

c 

ho  =  ho  /  0.005710147 
tshl  =  tshl  *  0.0254 
cshl  =  cshl  *  2.989067  /  12.0 
rhocs  =  rhocs  *  27679.9 
rcwt  =  rcwt  *  0.45359237 
forcig  =  forcig  *  2.989067 
covi  =  covi  /  27679.9 
chwi  =  chwi  *  0.45359237 
endif 

tmpi  =  0.0 

do  220  i  =  1,  nprop 

tmpi  =  tmpi  +  chwp(i) 
kpr  =  i 

call  prf 710 (pdp(i),  gdiap(i),  glenp(i),  nperfs(i),  0.,  frac(i) 
1  ,  volp(i),  surf(i),  dsdx(i)) 

tng(i)  =  chwp(i)  /  rhop(i)  /  volp(i) 
surfo(i)  =  surf(i) 
write(6,  1150)i,  tng(i) 

220  continue 

tmpi  =  tmpi  +  chwi 
write  (6,  1050 ) 
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read(2,  *,  end  =  790,  err  =  800)deltat,  deltap,  tstop,nzpi 

write  (6,  1190,  err  =  800)deltat,  deltap,  tstop 

write  (  *,  1200) 

deltat  =  deltat  *  0.001 

deltap  =  deltap  *  0.001 

tstop  =  tstop  *  .001 

if (igrad. eq. 2 .or . igrad.eq . 4 .or . igrad. eq. 6) go  to  230 
bore  =  (glr  *  grve  *  grve  +  aland  *  aland)  /  (glr  +  1 . ) 
bore  =  sqrt (bore) 

230  areab  =  pi  *  bore  *  bore  /  4. 
areaa=pi* (btdia/2 . ) **2 
areaba=areab-areaa 

lambda  =  1.  /  ({13.2  +  4.  *  loglOdOO.  *  bore))  **  2) 

iplot  =  0 

pltdt  =  deltat 

pltt  =  0. 

pmaxm  =  0.0 

pmaxbr  =  0.0 

pmaxba  =  0.0 

tpmaxm  =  0.0 

tpmxbr  =  0.0 

tpmxba  =  0.0 

tpmax  =  0.0 

a(l)  =  0.5 

a(2)  1 ,  -  sqrt(2.)  /  2. 

a(3)  =  1.  +  sqrt(2.)  /  2, 

a(4)  =  1.  /  6. 

b(l)  =  2. 

b(2)  =  1. 

b(3)  -  1. 

b(4)  =  2. 

ak(l)  =  0.5 

ak(2)  =  a(2) 

ak(3)  =  a(3) 

ak(4)  =  0.5 


240 


vpO  =  0  .  C 
trO  =  0.0 
tew  “  0.0 

if  (igrad.eq.  3)  chmlen  =  chaun  /  areab 
if (igrad.eq. 5) chmlen= (cham+btvol) /areab 
zb  =  chmlen 
zp  =  chmlen 


grlen  =  0. 
grdiam  =  0 . 
egama  =  0 . 
do  240  i  =  1,  nprop 
grlen  «  grlen  + 
grdiam  =  grdiam 
ibo ( i )  =0 
egama  =  egama  + 
nsl(i)  =  0 
vpO  =  chwp(i)  / 
continue 

volgi  =  cham  -  vpO  -  chwi  *  covi 
grlen  =  grlen  /  (tmpi  -  chwi) 
grdiam  =  grdiam  /  (tmpi  -  chwi) 
egama  =  (egama  +  chwi  *  gamai)  /  tmpi 
ism  =  0 


chwp(i)  * 
+  chwp(i) 

chwp ( i )  * 

rhop(i)  + 


glenp(i) 

*  gdiap(i) 

gamap ( i ) 

vpO 
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odlnr  =  0 . 

vfO  *  cham  -  vpO 

epsO  =  1 .  -  vpO  /  cham 

eps  =  epsO 

gasden  =  chwi  /  vfO 

prden  =  tmpi  /  vpO 

ug  =  0 . 

up  =  0 . 

pmean  =  forcig  *  chwi  /  volgi 
if (ihl . eq. 1 ) pmean  =  forcig  *  chwi  / 

1  (highv  -  vpO  -  chwi  *  covi) 
phase  «  pmean 
pbrch  =  pmean 
opbase  =  pmean 
volg  =  volgi 
volgi  =  volgi  +  vpO 
wallt  =  twal 
tgas  =  tempi 
told  =  0. 
tgaso  =  tgas 
dtgaso  =  0 . 
covl  =  covi 
t  =  0. 
ptime  =  0.0 
ibrp  =  12 
z(3)  =  1. 

nde  =  ibrp  +  nprop 

if (mode .eq. 1) write (6,  1210)areab,  pmean,  vpO,  volgi 
if (mode . eq . 2 ) then 

argl  =  areab  /  0.00064516 
arg2  =  pmean  /  6894.757 
arg3  =  vpO  /  0.000016387064 
arg4  =  volgi  /  0.000016387064 
write  (6,  1220) argl,  arg2,  arg3,  arg4 
endif 

write(6,  1236) title, vsn 
write(6,  1230) 
if (mode .eq. 1) write (6,  1232) 
if (mode .eq. 2) write (6,  1234) 
lines  =  4 
1 inmax  =  62 
iswl  =  0 
prwtO  =  prwt 
250  continue 

do  690  j  =  1,  4 
c 

c  find  barrel  resistance 

c 

if (ivpm.ne. 1) go  to  270 
do  260  k  =  2,  nvpmp 

if(y(2)  +  y  (7)  , It .pro jtr (k) ) go  to  270 
prwt  =  pro jms (k) 

260  continue 

270  if (ihl . eq. 1) go  to  300 

do  280  k  =  2,  npts 

if(y(2)  +  y (7) . ge . trav(k) ) go  to  280 
go  to  290 
280  continue 
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k  =  npts 

290  resp  =  (trav(k)  -  y(2)  -  y{7))  /  (trav(k)  -  trav(k  -  1) ) 

resp  =  br(k)  -  resp  *  (br(k)  -  br(k  -  1)) 

c 

c  find  mass  fraction  burned 

c 

300 


1 


310 


320 
c 

c  energy  loss  to  projectile  translation 
c 

elpt  =  y{ll) 
c 

c  elpt“prwt*y (1) *y  (1) /2 . 

c 

eptdot  =  prwt  *  y(l)  *  z(l) 
z(ll)  =  eptdot 
c 

c  energy  loss  due  to  projectile  rotation 
c 

elpr  =  y(12) 

elpr=pi*pi*prwt* ( (y (1) +y (6) ) **2) /4 . *twst*twst 
c 

eprdot  =  pi  *  pi  *  prwt  *  (y(l)  +  y(6))  *  (z(l)  +  z(6)) 
1  /  2 .  *  twst  *  twst 

z(12)  =  eprdot 
c 

c  energy  loss  due  to  gas  and  propellant  motion 

c 

if (igrad. eq. 1) go  to  340 
if {igrad.eq.3)go  to  350 
if (igrad. eq. 4 ) go  to  330 
if (igrad. eq. 5) go  to  352 
if (igrad. eq. 6)go  to  355 
pt  =  y(2)  +  y(7) 


do  320  k  =  1,  nprop 
kpr  =  k 

if (ibo (k) .eq. 1) goto 320 
nsll  =  0 

call  prf710 (pdp (k) ,  gdiap(k),  glenp(k),  nperfs(k), 
y(ibrp  +  k) ,  frac(k),  volp{k),  surf(k),  dsdx(k)) 
nsl(k)  =  nsll 
if (nsl (k) .eq. 0) goto  310 
if (nslp (k) .eq . 1 ) go  to  310 
write(6,  1240)k 
lines  =  lines  +  1 
nslp(k)  «=  1 
tsl  (k)  =  y (3) 
ism  =  1 
continue 

if (frac  (k) . It . . 9999)  go  to  320 

frac(k)  =  1. 

tbo(k)  =  y(3) 

ibo(k)  =  1 

ism  =  1 

write(6,  1250)k 
lines  =  lines  +  l 
continue 

if (ihl.eq.l)goto370 
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vzp  =  bvol  +  areab  *  pt 

j4zp  =  bint (4)  +  ((bvol  +  areab  *  pt)  **  3  -  bvol  **  3)  /  3. 

/  areab  /  areab 

elgpm  =  tmpi  *  y(l)  *  y(l)  *  areab  *  areab  *  j4zp  /  2.  /  vzp 
1  /  vzp  /  vzp 

go  to  360 

330  pb  =  y(7)  +  y(10) 

vzb  “  bvol  +  areab  *  pb 

j4zb  =  bint (4)  +  (vzb  **  3  -  bvol  **  3)  /  3.  /  areab  /  areab 
elgpm  =  (1.  -  eps)  *  up  *  up  *  areab  **  2  *  prden  *  j4zb  +  eps 
1  *  ug  *  ug  *  areab  **  2  *  gasden  *  j4zb 

elgpm  =  elgpm  /  2.  /  vzb  /  vzb  +  gasden  *  areab  *  ullen  /  6. 

1  *  (3.  *  y(l)  *  yd)  +  3.  *  yd)  *  ullen  *  dlnrho  +  ullen  **  2 

1  *  dlnrho  **  2) 

c 

c  approximate  epdot 
c 

epdot  =  tmpi  *  yd)  *  z(l)  /  3. 
go  to  360 

340  elgpm  =  tmpi  *  {y(l)  *  y(l)  -  y(l)  *  y(6)  +  y(6)  *  y(6))  /  6. 

go  to  360 

350  elgpm  =  areab  *  zb  /  6.  *  (eps  *  gasden  *  ug  *  ug  +  (1.  -  eps) 
1  *  prden  *  up  *  up) 

elgpm  =  elgpm  +  gasden  *  areab  *  ullen  /  6.  *  (3.  *  yd)  * 

1  y(l)  +  3.  *  y(l)  *  ullen  *  dlnrho  +  ullen  **  2  *  dlnrho  **  2) 
c 

approximate  epdot 

epdot  =  tmpi  *  y(l)  *  z(l)  /  3. 
go  to  360 

352  zp=y (2) +y (7 ) +chmlen 

za=zp-btlen 
vzp=zp*areab-btvol 
vza=za*areab 
c  write (6, 2050) tmpi, vzp, areaa, y (1 ) 

c  elgpm=tmpi*y (1) *y (1) /vzp/2 . 

c  elgpm=elgpm* (areab*areab/3 . /vzp/vzp* (areab*za**3 

c  &  + ( (areab* za+areaba* (zp-za) ) **3- (areab* za) **3) /areaba/ areaba) 

c  &  -areaa*areab/vzp/areaba/areaba* ( (areab*za+areaba* (zp-za) ) 

c  &  **2- (areab*za) **2) +areaa*areaa/areaba* ( zp-za) ) 

elgpml=tmpi*areab**2*y (l)**2/2./vzp/vzp/vzp 
taql= (areab*za**3/3.+ 

&  (areab*za+ (zp-za) *areaba) **3/3. /areaba/areaba- (areab*za) **3/ 

&  3 . /areaba/areaba) 

elgpm2=tmpi*areab*areaa*y (1) **2/vzp/vzp 
taq2= ( (areab*za+ (zp-za) 

&  *areaba) **2/2 . /areaba/areaba- (areab*za) **2/2 . /areaba/areaba) 

elgpm3=tmpi*y (1 ) **2*areaa**2/2 . /vzp 
taq3= (zp-za) /areaba 

elgpm=elgpml*taql-elgpm2*taq2+elgpm3*taq3 
c  write (6, 2050) tmpi, vzp, areaa, y (1) 

c2050  format('  tmpi' , el7 . 10, '  vzp' , el7 . 10 , '  areaa' , el7 . 10, '  y(l)' 

c  &  , el7 . 10 ) 

go  to  360 
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355  pt=y(2)+y(7) 

ptl=pt+chmlen 

call  jint (btdia,btlen,ptl,ptl,nchpts,chdist,chdiam,bint,bvolzp) 

vzp=bvolzp 
qlzp=bint (10) 
q2zp=bint (2 ) 
q6zp=bint  (4 ) 
q8zp=bint  (8 ) 
q3zp=bint (5) 
q4zp=bint (6) 
q5zp=bint (7 ) 
q9zp=bint (9) 

write  (6,  ■76)qlzp,q2zp,q3zp,q4zp,q5zp,q62p,q7zp,q8zp,ptl,pt2 

76  format  (lx,  lOell  .4)  r-  r-  r 

delta=l . 

ptl=chmlen+y (2) +y (7) 
pt2=chmlen+y (2) +y (7) -btlen 
zp=ptl 
za=pt2 

call jint (btdia, btlen, ptl, pt2, nchpts, chdist, chdiam, bint, bvolza) 

vza  =  bvolza 

qlza=bint (1 ) 

q22a=bint (2) 

q4za=bint (6) 

q5za=bint  (7) 

q6za=bint (4 ) 

q7za=bint (3) 

q3zpza=q3zp 

q9zpza=qlzp 

c  write(6,77)qlza,q2za,q3za,q4za,q5za,q6za,q7za,q8za,ptl,pt2 

77  rormat  (lx,  lOell  .4)  ft'  ft' 

c  elgpm=tmpi*y(l)*y(i)*areab**2*q6zp/2./vzp**3  -  tmpi*areaa*areab* 

C  &y(l)*y(l)*q9zpza/vzp/vzp  +  tmpi*y (l) *y (l) *areaa**2*q3zpza/ 

c  &2,/vzp 

elgpml=tmpi*y (1) *y(l) *areab**2*q6zp/2 . /vzp**3 
elgpm2=tmpi*areaa*areab*y (l)*y{l) *qlzp/vzp/vzp 
elgpm3-tmpi*y (1) *y (1) *areaa**2*q3zp/2 . /vzp 
elgpm=elgpml-elgpm2+elgpm3 
go  to  360 


energy  loss  from  bore  resistance 

360  elbr  =  y(4) 

z(4)  =  areab  *  resp  *  (y(l)  +  y(6)) 
ebrdot  =  z(4) 

c  energy  loss  due  to  recoil 
c 

elrc  =  rcwt  *  y(6)  *  y(6)  /  2. 
erdot  =  rcwt  *  y(6)  *  z(6) 
c 

c  energy  loss  due  to  heat  loss 
c 

areaw  =  cham  /  areab  *  pi  *  bore  +  2 
1  (y(2)  +  y(7)) 

370  avden  =  0.0 
avc  =  0.0 


*  areab  +  pi  *  bore  * 


/  (gcimap(k) 


zl9 

f rac (k) 


c 

c 

c 

c 


c 

c 

c 


avcp  =  0.0 
zl8  =  0 
zl9  =  0 

do  380  k  =  1,  nprop 

zl8  =  forcp(k)  *  gamap(k)  *  chwp(k)  *  frac{k) 

1  -  1.)  /  tempp(k)  +  zl8 

zl9  =  chwp(k)  *  frac(k)  + 
avden  =  avden  +  chwp(k)  * 

380  continue 

avcp  =  (zl8  +  forcig  *  gamai  *  chwi  /  (gamai  -  1.)  /  tempi)  / 
1  {zl9  +  chwi) 

avden  =  (avden  +  chwi)  /  (volg  +  covl) 

avvel  =  .5  *  (y(l)  +  y{6)) 

htns  =  lambda  *  avcp  *  avden  *  avvel  +  ho 

z{5)  =  areaw  *  htns  *  (tgas  -  wallt)  *  hi 

elht  =  y(5) 

ehdot  =  z (5) 

wallt  =  (elht  +  htfr  *  elbr)  /  cshl  /  rhocs  /  areaw  /  tshl  + 

1  twal 

energy  loss  due  to  air  resistance  in  tube 

(assume  no  drag  resistance  on  air/tube  interface) 

if (ihl . eg. 1) goto  410 
air  =  iair 

z(8)  =  y(l)  *  pgas  *  air 
elar  =  areab 
eddot  =  z (8 ) 

recoil 


y(8) 

*  areab 


z(6)  =  0.0 

if (pbrch . le . rp ( 1)  /  areab)go  to  400 
rfor  =  rp(2) 

if(y(3)  -  trO . ge . tr (2) ) go  to  390 


390 


400 

410 


-  trO) )  /  (tr  (2)  -  tr  (1) ) 
(rp(2)  -  rp(l)  ) 

(pbrch  -  rfor  /  areab  -  resp) 
0.0 


c 

c 

c 


rfor  =  (tr(2)  -  (y(3) 
rfor  =  rp(2)  -  rfor  * 
z(6)  =  areab  /  rcwt  * 
if (y (6)  .  It .  0 . 0) y (6)  = 
z  (7)  =  y  (6) 
goto  410 
trO  =  y(3) 
continue 

calculate  gas  temperature 

eprop  =  0.0 
rprop  =  0.0 
drnfogt  =  0.0 
dmfog  =  0.0 
do  420  k  =  1,  nprop 

eprop  =  eprop  +  forcp(k)  *  chwp(k)  *  frac(k)  /  (gamap(k) 

L  -  1.) 

rprop  =  rprop  +  forcp(k)  *  chwp(k)  *  frac(k)  /  (gamap(k) 

L  -  1.)  /  tempp(k) 

drnfogt  =  drnfogt  +  forcp(k)  *  rhop(k)  *  tng(k)  *  surf(k)  * 
L  z(ibrp  +  k)  /  ((gamap{)c)  -  1.)  *  tempp(k)) 

dmfog  =  dmfog  +  forcp(k)  *  rhop(k)  *  tng(k)  *  surf(k)  * 
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1  z(ibrp  +  k.)  /  (gamap<k)  -  1.) 

420  continue 

tenerg  =  elpt  +  elpr  +  elgpm  +  elbr  +  elrc  +  elht  +  elar 

tgas  =  (eprop  +  forcig  *  chwi  /  (gamai  -  1.)  -  elpt  -  elpr  - 

1  elgpm  -  elbr  -  elrc  -  elht  -  elar)  /  (rprop  +  forcig  *  chwi 

1  /  (gamai  -  1.)  /  tempi) 

tedot  =  epdot  +  eprdot  +  eddot  +  ebrdot  +  erdot  +  ehdot+eptdot 
dtgas  =  (dmfog  -  tedot  -  tgas  *  dmfogt)  /  (rprop  +  forcig  * 

1  chwi  /  (gamai  -  1.)  /  tempi) 

c 

c  find  free  volume 

c 

vl  =  0.0 

covl  =  0.0 

do  430  k  =  1,  nprop 

vl  =  chwp(k)  *  (1.  -  frac(k))  /  rhop(k)  +  vl 
covl  =  covl  +  chwp(k)  *  covp(k)  *  frac(k) 

430  continue 

volg  =  volgi  +  areab  *  (y(2)  +  y(7))  -  vl  -  covl 
if ( ihl . eg . 1 ) volg  =  highv  -  vl  -  covl 
c 

c  calculate  mean  pressure 

c 

rl  =  0.0 

do  440  k  =  1,  nprop 

rl  =  rl  +  forcp(k)  *  chwp(k)  *  frac(k)  /  tempp(k) 

440  continue 

pmean  =  tgas  /  volg  *  (rl  +  forcig  *  chwi  /  tempi) 
if ( ihl . eq . 1 ) go  to  640 
resp  =  resp  +  pgas  *  air 
if ( igrad. eq. 1 ) go  to  590 
if (igrad.eq.2)go  to  450 
if (igrad. eq.3)go  to  470 
if ( igrad. eq . 4 ) go  to  540 
if ( igrad. eq . 5 ) go  to  582 
if ( igrad . eq . 6) go  to  585 
450  if ( iswl . ne . 0 ) go  to  460 

phase  =  pmean 
pbrch  =  pmean 

if (phase . gt . resp  +  l.)iswl  =  1 
go  to  620 
c 

c  use  chambrage  pressure  gradient  equation 

c 

460  jlzp  =  bint(l)  +  (bvol  *  pt  +  areab  /  2.  *  pt  *  pt)  /  areab 

j2zp  =  (bvol  +  areab  *  pt)  **  2  /  areab  /  areab 
j3zp  =  bint (3)  +  areab  *  bint(l)  *  pt  +  bvol  *  pt  *  pt  /  2 .  + 
1  areab  *  pt  *  pt  *  pt  /  6 . 

a2t  =  -  tmpi  *  areab  *  areab  /  prwt  /  vzp  /  vzp 
alf  =  1.  -  a2t  *  jlzp 

alt  =  tmpi  *  areab  *  (areab  *  y(l)  *  y(l)  /  vzp  +  areab  *  resp 
1  /  prwt)  /  vzp  /  vzp 

bt  =  -  tmpi  *  y(l)  *  y(l)  *  areab  *  areab  12.1  vzp  /  vzp/vzp 
bata  =  -  alt  *  jlzp  -  bt  *  j2zp 
gamma  =  alf  +  a2t  *  j3zp  /  vzp 
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c 

c 

c 

c 

c 

c 


c 


delta  =  bata  +  alt  *  j3zp  /  vzp  +  bt  *  j4zp  /  vzp 

calculate  base  pressure 

phase  =  (pmean  -  delta)  /  geunma 

calculate  breech  pressure 

pbrch  =  alf  *  phase  +  bata 
go  to  610 

use  2  phase  gradient  equation 


470  if (iswl . ne . 0) goto  480 
phase  =  pmean 
pbrch  =  pmean 

if (phase . gt . resp  +  l)iswl  =  1 
go  to  620 

480  if (iswl .eq. 2) go  to  580 

vzp  =  Cham  +  areab  *  (y(2)  +  y(7)) 
vzb  =  Cham  +  areab  *  (y(10)  +  y(7)) 
phi  =  C 
phidot  =  0 . 
dmorho  =  0 . 
dmcov  =  0 . 
dmromw  =  0 . 
rmomw  =  0 . 
vfree  =  vzp  -  vl 
do  490  k  =  1,  nprop 

rmomw  =  rmomw  +  chwp(k)  *  frac(k)  *  forcp(k)  /  tempp(k) 
phi  =  chwp(k)  *  frac(k)  +  phi 
if (ibo (k) .eq. 1) go  to  490 

dmorho  =  dmorho  +  tng(k)  *  surf(k)  *  z(ibrp  +  k) 
phidot  =  rhop(k)  *  tng(k)  *  surf(k)  *  z(ibrp  +  k)  +  phidot 
dmcov  =  rhop(k)  *  tng(k)  *  surf (k)  *  z(ibrp  +  k)  *  covp(k) 
1  +  dmcov 

dmromw  =  dmromw  +  rhop(k)  *  tng(k)  *  surf  (k)  *  z(ibrp  +  k) 
1  *  forcp(k)  /  tempp(k) 

490  continue 

rmomw  =  rmomw  +  chwi  *  forcig  /  tempi 
gasmas  =  phi  +  chwi 
gasden  =  gasmas  /  vfree 
phi  =  (phi  +  chwi)  /  tmpi 
if  (phi . gt . 0 . 999)  then 
iswl  =  2 

rbm  =  phase  /  pmean 
rbrm  =  pbrch  /  pmean 
if (phi . ge . 1 . ) go  to  580 
endif 

dmdt  =  phidot 
phidot  =  phidot  /  tmpi 

vdotov  =  (dmorho  +  areab  *  y(l))  /  vfree 
dlnrho  =  dmdt  /  gasmas  -  vdotov 
dvoldt  =  dmorho  +  areab  *  y(l)  -  dmcov 
c 

c  get  time  derivative  of  mean  pressure 

c 

dpmdt  =  (dmromw  *  tgas  -  pmean  *  dvoldt  +  dtgas  *  rmomw)  /volg 
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volprp  *  0. 
effdia  =  0. 
dmdmdt  =  0 . 
dmdinor  =  0 . 
avelen  =  0. 
avedia  =  0. 
do  500  k  =  1,  nprop 

if (ibo (k) .eq. 1) go  to  500 

volprp  =  volprp  +  (1.  -  frac(k))  *  chwp(k)  /  rhop(k) 
dmdmdt  “  dmdmdt  +  rhop(k)  *  tng(k)  *  dsdx(k)  *  z(ibrp  +  k) 
1  *  z (ibrp  +  k) 

dmdmdt  -  dmdmdt  +  rhop(k)  *  tng(k)  *  surf(k)  *  d2xdt2(k) 
dmdmor  =  dmdmor  +  (dsdx(k)  *  z(ibrp  +  k)  •  *  2  surf(k)  * 
1  d2xdt2(k))  *  tng(k) 

effdia  =  effdia  +  6.  *  volp(k)  /  surf  (k)  *  (1.  -  frac(k)) 
1  *  chwp(k) 

500  continue 

clt  =  dmdmdt  /  gasmas  -  dn  jmor  /  vfree  +  vdotov  **  2  -  (dmdt 
1  /  gasmas)  **  2 

d21nr  =  clt  -  areab  **  2  *  phase  /  vfree  /  prwt 

d21nr  =  d21nr  +  areab  *  areab  *  resp  /  vfree  /  prwt 

zp  =  chmlen  +  y(2)  +  y(7) 

zb  =  chmlen  +  y(10)  +  y{7) 

ullen  =  zp  -  zb 

enow  =  tmpi  -  gasmas 

vp  =  y(l) 

effdia  =  effdia  /  enow 
prden  =  enow  /  volprp 
up  =  y(9) 

phistr  =  phi  -  gasden  *  areab  *  ullen  /  tmpi 
ulldot  »  vp  -  up 

dphist  =  phidot  -  gasden  *  areab  /  tmpi  *  (ulldot  +  ullen  * 

1  dlnrho) 

eps  =1.  -  (1.  -  phi)  ♦  tmpi  /  prden  /  vzb 

epsdot  =  phidot  *  tmpi  /  prden  /  vzb  +  (1.  -  phi)  *  tmpi  *  up 
1  *  areab  /  prden  /  vzb  /  vzb 

ug  =  up  +  (vp  +  ullen  *  dlnrho  -  up)  /  eps 
alam  =  (1.5  *  grlen  /  grdiam)  **  .666666667 
alam  =  (0.5  +  grlen  /  grdiam)  /  alam 
alam  =  alam  **  2.17 
c 

c  vis  kg/s/m 
c 

vis  =  .00007 

ren  =  gasden  /  vis  *  effdia  *  abs (ug  ~  up) 
if (ren . It . 1 . ) ren  =  1. 

fsrg  =  2.5  *  alam  /  ren  **  .081  *  ((1.  -  eps)  /  (1.  -  epsO)  * 
1  epsO  /  eps)  **  .45 
fsc  =  fsrg  *  fsO 

phi2  =  1,  -  phi  -  phistr  *  (1.  -  eps)  /  eps 

philp  =  dphist  *  ug  -  phidot  *  up  -  phistr  *  epsdot  /  eps  / 

1  eps  *  (vp  +  ullen  *  dlnrho  -  up)  +  phistr  *  ulldot  *  dlnrho 

1  /  eps  +  2.  *  phistr  *  ug  /  zb  *  (ug  -  up) 

philp  =  philp  +  phi2  *  gasden  /  effdia  /  prden  * 

1  (ug  -  up)  **  2  *  fsc 

ak2  =1.  /  (1.  -  phi2  *  tmpi  /  prden  /  vzb) 
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c 

c 


phil  =  philp  +  phistr  *  z(l)  /  eps  +  ullen  *  phistr  *  d21nr  / 
1  eps 

acceleration  of  forward  boundary  of  propellant  bed 

z(9)  =  gasden  *  (ug  -  up)  **  2  *  fsc  /  prden  /  effdia  +  tmpi 
1  *  phil  *  ak2  /  vzb  /  prden 

z(10)  =  y(9) 

e  =  phistr  /  eps  *  (1.  -  ullen  *  areab  /  vfree)  *  areab  /  prwt 
dd  =  ullen  *  phistr  *  clt  /  eps 
akll  =  tmpi  *  e  *  ak2  /  zb  /  vzb 

akl2  =  tmpi  *  a>.2  *  (philp  +  dd)  /  zb  /  vzb  -  akll  *  resp 
phase  =  pmean  -  akl2  *  zb  *  zb  /  2 .  +  gasden  *  ullen  *  resp  * 
1  areab  /  prwt 

phase  =  phase  +  akl2  *  zb  *  zb  *  (zb  /  3.  +  ullen)  /  2.  /  zp 

phase  =  phase  -  gasden  *  ullen  **  2  *  areab  *  resp  /  2.  /  zp 

1  /  prwt 

phase  =  phase  -  gasden  *  ullen  **2/2.  *  (1.  -2.  *  ullen  / 
1  3.  /zp)*(clt-  dlnrho  **  2) 

phase  =  phase  -  areab  **  2  *  gasden  *  ullen  **  2  *  (1.-2.  * 
1  ullen  /  3.  /  zp)  *  resp  /  prwt  /  vfree  /  2. 

deno  =  -  akll  *  zb  **  3  /  6.  /  zp  -  ullen  *  akll  *  zb  *  zb  / 

1  2 .  /  zp 

deno  =  deno  +  gasden  *  ullen  *  areab  /  prwt  -  areab  **  2  * 

1  gasden  *  ullen  **  2  *  (1.  -  2.  *  ullen  /  3.  /  zp)  /  2.  / 

1  vfree  /  prwt 

deno  =  deno  -  gasden  *  ullen  **  2  *  areab  12.  /  zp  /  prwt  + 

1  1.  +  akll  *  zb  *  zb  /  2. 

phase  =  phase  /  deno 
if (ism.eq.0)goto530 
if ( ism. eq . 1 ) gotoSlO 
goto520 
510  ism  =  2 

tss  =  sqrt(egama  *  rmomw  /  gasmas  *  tgas) 
write(6,  *  )tss 

tss  =  ullen  /  (ullen  *  odlnr  +  tss) 
tso  =  y(3) 

write(6,  *  )tss,  tso 

520  coefbp  =  (tss  +  tso  -  y(3)  -  deltat)  /  tss 

if (coefbp . gt . 1 . ) coefbp  =  1. 
if (coefbp . le . 0 . ) then 
coefbp  =  0. 
ism  =  0 
endif 

phase  =  coefbp  *  opbase  +  (1.  -  coefbp)  *  phase 
if (mode .eq. 1) write (6,  *  ) coefbp,  opbase,  phase,  ism 
if (mode . eq . 2 ) then 

argl  =  opbase  /  6894.757 
arg2  =  phase  /  6894.757 
write  (6,  *  ) coefbp,  argl,  arg2,  ism 
endif 

530  odlnr  =  dlnrho 

opbase  =  phase 

pbrch  =  phase  *  (1.  +  akll  *  zb  *  zb  /  2 .  +  gasden  *  ullen  * 

1  areab  /  prwt  -  areab  **  2  *  gasden  *  ullen  **  2  /  2.  /  vfree 

1  /  prwt) 
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pbrch  =  pbrch  +  akl2  *  zb  *  zb  /  2.  -  gasden  *  ullen  *  areab 
1  *  resp  /  prwt 

pbrch  =  pbrch  +  gasden  *  ullen  **2/2.  *  (clt  -  dlnrho  **  2) 

pbrch  =  pbrch  +  areab  **  2  *  gasden  *  ullen  **  2  *  resp  /  2. 

1  /  vfree  /  prwt 

go  to  610 
c 

c  using  rga  gradient 
c 

540  if (iswl . ne . 0) go  to  550 
phase  =  pmean 
pbrch  =  pmean 

if (phase. gt .resp  +  l.)iswl  =  1 
go  to  620 

550  if (iswl .eq. 2) go  to  580 

vzp  =  Cham  +  areab  *  {y(2)  +  y(7)) 
vzb  =  Cham  +  areab  *  (y{10)  +  y(7)) 

jlzb  =  bint(l)  +  (bvol  *  pb  +  areab  /  2.  *  pb  *  pb)  /  areab 

j2zb  =  (bvol  +  areab  *  pb)  **  2  /  areab  /  areab 

j3zb  =  bint (3)  +  areab  *  bint(l)  *  pb  +  bvol  *  pb  *  pb  /  2 .  + 

1  areab  /  6.  *  pb  **  3 

phi  =  0. 
phidot  =  0. 
dmorho  =  0. 
dmcov  =  0 . 
dmromw  =  0 . 
rmomw  =  0 . 
vfree  =  vzp  -  vl 
do  560  k  =  1,  nprop 

rmomw  =  rmomw  +  chwp(k)  *  frac(k)  *  forcp(k)  /  tempp(k) 
phi  =  chwp(k)  *  frac(k)  +  phi 
if (ibo (k) .eq. 1) go  to  560 

dmorho  =  dmorho  +  tng(k)  *  surf{k)  *  z(ibrp  +  k) 
phidot  =  rhop(k)  *  tng(k)  *  surf(k)  *  z(ibrp  +  k)  +  phidot 
dmcov  =  rhop(k)  *  tng(k)  *  surf(k)  *  z(ibrp  +  k)  *  covp(k) 
1  +  dmcov 

dmromw  =  dmromw  +  rhop(k)  *  tng(k)  *  surf  (k)  *  z(ibrp  +  k) 
1  *  forcp(k)  /  tempp(k) 

560  continue 

rmomw  =  rmomw  +  chwi  *  forcig  /  tempi 
gasmas  =  phi  +  chwi 
gasden  =  gasmas  /  vfree 
phi  =  (phi  +  chwi)  /  tmpi 
if  (phi .gt . 0 . 99)  then 
iswl  =  2 

rbm  =  phase  /  pmean 
rbrm  =  pbrch  /  pmean 
if (phi . ge . 1 . ) go  to  580 
endif 

dmdt  =  phidot 
phidot  =  phidot  /  tmpi 

vdotov  =  (dmorho  +  areab  *  y(l))  /  vfree 
dlnrho  -  dmdt  /  gasmas  -  vdotov 
dvoldt  =  dmorho  +  areab  *  y(l)  -  dmcov 
c 

c  get  time  derivative  of  mean  pressure 
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dpmdt  =  (dmromw  *  tgas  —  pmean  *  dvoldt  +  dtgas  *  rmomw)  /volg 

volprp  =  0 . 

effdia  =  0. 

dmdmdt  =  0 . 

dmdmor  =  0 . 

avelen  =  0. 

avedia  =  0 . 

do  570  k  =  1,  nprop 

if (ibo (k) .eq. 1) go  to  570 

volprp  =  volprp  +  (1.  -  frac(k))  *  chwp(k)  /  rhop{k) 
dmdmdt  =  dmdmdt  +  rhop(k)  *  tng(k)  *  dsdx(k)  *  z(ibrp  +  k) 
1  *  z(ibrp  +  k) 

dmdmdt  =  dmdmdt  +  rhop(k)  *  tng{k)  *  surf  (k)  *  d2xdt2 (k) 
dmdmor  =  dmdmor  +  (dsdx(k)  *  z(ibrp  +  k)  **  2  +  surf (k)  * 
1  d2xdt2(k))  *  tng(k) 

effdia  =  effdia  +  6.  *  volp(k)  /  surf(k)  *  (1.  -  frac(k)) 
1  *  chwp(k) 

570  continue 

clt  =  dmdmdt  /  gasmas  -  dmdmor  /  vfree  +  vdotov  **  2  -  (dmdt 
1  /  gasmas)  **  2 

d21nr  =  clt  -  areab  **  2  *  phase  /  vfree  /  prwt 

d21nr  =  d21nr  +  areab  *  areab  *  resp  /  vfree  /  prwt 

zp  =  chmlen  +  y(2)  +  y(7) 

zb  =  chmlen  +  y(10)  +  y(7) 

ullen  =  zp  -  zb 

enow  =  tmpi  -  gasmas 

vp  =  y ( 1 ) 

effdia  =  effdia  /  enow 
prden  ■=  enow  /  volprp 
up  =  y  /  9 ) 

phistr  =  phi  -  gasden  *  areab  *  ullen  /  tmpi 
ulldot  =  vp  -  up 

dphist  =  phidot  -  gasden  *  areab  /  tmpi  *  (ulldot  +  ullen  * 

1  dlnrho) 

eps  =1.  -  (1.  -  phi)  *  tmpi  /  prden  /  vzb 

epsdot  =  phidot  *  tmpi  /  prden  /  vzb  +  (1.  -  phi)  *  tmpi  *  up 
1  *  areab  /  prden  /  vzb  /  vzb 

ug  =  up  +  (vp  +  ullen  *  dlnrho  -  up)  /  eps 
alam  =  (1.5  *  grlen  /  grdiam)  **  .666666667 
alam  =  (0.5  +  grlen  /  grdiam)  /  alam 
alam  =  alam  **  2.17 
c 

c  vis  kg/s/m 

c 

vis  =  .00007 

ren  =  gasden  /  vis  *  effdia  *  abs (ug  -  up) 
if (ren. It . 1 . ) ren  =  1. 

fsrg  =  2.5  *  alam  /  ren  **  .081  *  ((1.  -  eps)  /  (1.  -  epsO)  * 
1  epsO  /  eps)  **  .45 
fsc  =  fsrg  *  fsO 

phi2  =  1.  -  phi  -  phistr  *  (1.  -  eps)  /  eps 

philp  =  dphist  *  ug  -  phidot  *  up  -  phistr  *  epsdot  /  eps  / 

1  eps  *  (vp  +  ullen  *  dlnrho  -  up)  +  phistr  *  ulldot  *  dlnrho 

1  /  eps  +2.  *  areab  *  phistr  *  ug  /  vzb  *  (ug  -  up) 

philp  =  philp  +  phi2  *  gasden  /  effdia  /  prden  * 
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(ug  -  up)  **2  *  fsc 

ak2  =1.  /  (1.  -  phi2  *  tmpi  /  prden  /  vzb) 

phil  =  philp  +  phistr  *  z(l)  /  eps  +  ullen  *  phistr  *  d21nr  / 
1  eps 

c 

c  acceleration  of  forward  boundary  of  propellant  bed 

c 

z(9)  =  gasden  *  (ug  -  up)  **  2  *  fsc  /  prden  /  effdia  +  tmpi 
1  *  phil  *  a)c2  /  vzb  /  prden 

z(10)  =  y(9) 

phi3  =  phistr  *  ug  *  ug  +  (1.  -  phi)  *  up  *  up 
e  =  1 .  -  ullen  *  areab  /  vf ree 
dd  =  ullen  *  phistr  *  clt  /  eps 

alt  =  tmpi  *  areab  /  vzb  /  vzb  *  (phi3  *  areab  /  vzb  -  (philp 
1  +  dd  -  e  *  phistr  *  areab  *  resp  /  eps  /  prwt)  *  ak2) 

a2t  =  (  -  tmpi  *  e  *  phistr  *  areab  **  2  /  vzb  /  vzb  /  eps  / 

1  prwt)  *  ak2 

bt  =  -  tmpi  *  phi3  *  areab  **2/2.  /  vzb  /  vzb  /  vzb 
phase  =  pmean  -  gasden  *  ullen  **2/2.  *  (clt-  dlnrho  **  2 
1  +  areab  **  2  *  resp  /  prwt  /  vfree)  *  (1.  -  2.  *  areab  * 

1  ullen  /  3.  /  vzp) 

phase  =  phase  -  alt  *  j3zb  /  vzp  -  bt  *  j4zb  /  vzp  -  areab  * 

1  ullen  *  alt  *  jlzb  /  vzp 

phase  =  phase  -  areab  *  bt  *  ullen  *  j2zb  /  vzp  -  gasden  * 

1  areab  **  2  *  ullen  **  2  *  resp  /  2.  /  vzp  /  prwt  +  alt  * 

1  jlzb  +  bt  *  j2zb  +  areab  *  gasden  *  ullen  *  resp  /  prwt 

deno  =  1 .  areab  *  ullen  *  a2t  *  jlzb  /  vzp  -  gasden  *  areab 
1  **  2  *  ullen  **2/2.  /  vzp  /  prwt  +  a2t  *  j3zb  /  vzp  -  a2t 

1  *  jlzb  +  gasden  *  ullen  *  areab  /  prwt 

deno  =  deno  -  gasden  *  ullen  **  2  *  areab  **2/2.  /  vfree  / 
1  prwt  +  gasden  *  areab  **  3  *  ullen  **3/3.  /  vzp  /  vfree  / 
1  prwt 

phase  =  phase  /  deno 

pbrch  =  phase  *  (1.  -  a2t  *  jlzb  +  gasden  *  ullen  *  areab  / 

1  prwt  -  gasden  *  ullen  **  2  *  areab  **2/2.  /  vfree  /  prwt) 

1  +  gasden  *  ullen  **  2  /  2 .  *  (clt  -  dlnrho  **  2  +  areab  **  2 

1  *  resp  /  prwt  /  vfree)  -  alt  *  jlzb  -  bt  *  j2zb  -  areab  * 

1  gasden  *  ullen  *  resp  /  prwt 

go  to  610 

580  phase  =  rbm  *  pmean 

pbrch  =  rbrm  *  pmean 
go  to  610 

582  areazm=areab 

areaza=areazm-  pi  *  btdia**2/4. 
c  write ( G ,*) areazm,  areaza 

delta=l . 

qlzp= (vzp* *2- ( areab* za) **2) /2 . /areaba**2 

q2zp=vzp**2/areaba/areaba 

q3zp= (zp-za) /areaba 

q4zp=vzp/ areaba/areaba 

q5zp=l . /areaba/areaba 

q6zp= (areab*za**3/3 .+ (areab*za+areaba* (zp-za) ) **3/3 . /areaba 
&**2- ( areab* za) **3/3 . / areaba* *2) 
q8zp= (zp-za) **2/2 . 

q9zp= (vzp**3- ( areab *za) **3) /6 . /areaba* *2 
&  - (areab*za) **2* (zp-za) /2 ./areaba 
qlza=za**2/2 . 
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q2za=vza**2/areab/areab 
q7za=areab*za**3/ 6 . 
vp-y ( 1 ) 

bt=  -tmpi  *  y (1) *y (1) *areab*areab/2 . /vzp/vzp/vzp 

hl=  tmpi*areab*areaa*y (1) *y (1) / (vzp*vzp) 

akl=-tmpi*areab*vp*vp*vza*areaa/vzp/vzp/areaza**2 

akl=akl+tmpi*areaa**2*vp**2/vzp/areaza**2/2 . 

akl=akl+tmpi*vp*vp*areab**2*vza**2/2 . /vzp**3/areaza**2 

a jl=-tmpi*areaa**2*vp*vp/2 . /vzp 

fk=  -2 ,  *tnipi*vp**2*areazm*areaa/areaza/vzp* 

&  (areab*vza/vzp/areaztn-l . )  **2 
fk“  fk/ (areaza+areazm) 

rrl=  1.  +  tinpi*areab*areaa*qlza/prwt/vzp**2 
tal=tmpi*areab**2*y (1 ) **2/vzp**3 
ta4=tmpi*areab**2*resp/prwt/vzp**2 
zal=  (bt*q2za  +  (tal+ta4) *qlza) /rrl 
zaO=  l./rrl 

za2=  -  (tmpi*areab*areaba*qlza/prwt/vzp*’^2  ) /rrl 
a3  =  tmpi*areab**2*y (1) **2/vzp**3  -  tmpi*areab*areaa*zal/prwt/ 
&  vzp**2  +  tmpi*areab**2*resp/vzp**2/prwt 

a41=  -tmpi*areab*areaa*za2/vzp**2/prwt 
a42=  -tmpi*areab*areaba/vzp**2/prwt 
a4  =  a41+a42 

c  a4  =  -tmpi*areab*areaa*za2/vzp**2/prwt  -  tmpi*areab*areaba/ 

c  &  vzp**2/prwt 

a5  =  -  tmpi*areab*areaa*za0/prwt/vzp**2 
c3  =  tmpi*areaa**2*zal/prwt/vzp  -  tmpi*areaa*areab* 

&  resp/prwt/vzp  -  tmpi*areaa*areab*y (1) *y { 1 ) /vzp/vzp 
c41  =  tmpi*areaa*areaa*za2/vzp/prwt 
c42  =  tmpi*areaa*areaba/vzp/prwt 
c4  =  c41+c42 

c4  =  delta*tmpi*areaa*areaa*za2/vzp/prwt  +  cielta*tmpi*areaa* 

&  areaba/vzp/prwt 
c5  =  tmpi*areaa*areaa*zaO/vzp/prwt 

11  =  zal+fk+a3*qlzp+c3*q3zp+bt*q2zp+hl*q4zp+a jl*q5zp+akl 

12  =  l-2a2-a4 *qlzp  -  c4*q3zp 

13  =  zaO  +  a5*qlzp  +c5*q3zp 

bl  =  (a3*q7za+a3*q9zp+bt*q6zp+zal* (vzp-vza) +fk* (vzp-vza) 

*  +c3*q8zp+hl*qlzp+a jl*q3zp+akl* (vzp-vza) ) /vzp 
b2  =  (a4*q7za+a4*q9zp+za2* (vzp-vza) +c4*q8zp) /vzp 
b3  =  (vza+a5*q7za+a5*q9zp+za0* (vzp-vza) +c5*q8zp) /vzp 
c  calculate  base  pressure 

pbase=(pmean/b3  -  bl/b3  +  11/13) /(  12/13  +  b2/b3) 
c  calculate  breech  pressure 

pbrch=pmean/b3  -  bl/b3  -  pbase*b2/b3 
pza=zal  +  za2*pbase  +  zaO*pbrch 
c  write (6,7)y(3),z(l),y(l),y(2), pmean , phase , pbr ch 

c  write(6,7)  bl , b2, b3, 11, 12, 13 

c  write(6,7)  a3, a4 , a5, c3, c4 , c5 

z (1) = (areaa*pza+areaba*pbase-areab*resp) /prwt 
go  to  615 

585  If (za . gt . chdist (nchpts) ) goto  275 
Do  269  1=2, nchpts 

If (za . It . chdist (I) . and. za . gt . (chdist (I-l) ) ) goto  274 

269  continue 

274  diam  =(  za  -  chdist (I-l))  /  (  chdist (I)  -  chdist (I-l))  * 

&  (chdiam(I)  -  chdiam(I-l))  +  chdiam(I-l) 

areazm  =  pi  *  diam**2/4. 
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areaza  =•  areazm  -  pi  *  btdia**2/4. 

275  continue 
vp-y (1) 

bt-  -tmpi  *  y (1) *y (1) *areab*areab/2 . /vzp/vzp/vzp 

hl=  tmpi*areab*areaa*y (1) *y (1) / {vzp*vzp) 

akl=-tmpi*areab*vp*vp*vza*areaa/vzp/vzp/areaza**2 

akl=akl+tmpi*areaa**2*vp*  *2/vzp/areaza**2/2 . 

akl=akl+tmpi*vp*vp*areab**2*vza**2/2 , /vzp**3/areaza**2 

a jl=-tmpi*areaa**2*vp*vp/2 . /vzp 

fk=  -2 . *tmpi*vp**2*areazm*areaa/areaza/vzp* 

& (areab*vza/vzp/areazm-l . ) **2 
f k-  f k/ (areaza+areazm) 

rrl=  1.  +  tmpi*areab*areaa*qlza/prwt/vzp**2 
tal=tmpi*areab**2*y (1) **2/vzp**3 
ta4=tmpi*areab**2*resp/prwt /vzp**2 
zal=  (bt*q2za  +  (tal+ta4 ) *qlza) /rrl 
zaO=  l./rrl 

za2=  - (tmpi*areab*areaba*qlza/prwt/vzp**2  ) /rrl 
a3  =  tmpi*areab**2*y (1) **2/vzp**3  -  tmpi*areab*areaa*zal/prwt/ 
&  vzp**2  +  tmpi*areab**2*resp/vzp**2/prwt 

a41=  -tmpi*areab*areaa*za2/vzp^*2/prwt 
a42=  -tmpi*areab*aLeaba/vzp**2/prwt 
a4  =  a41+a42 

c  a4  =  -tmpi*areab*areaa*za2/vzp**2/prwt  -  tmpi*areab*areaba/ 
c  &  vzp**2/prwt 

a5  =  -  tmpi*areab*areaa*za0/prwt/vzp**2 
c3  =  tmpi*areaa**2*zal/prwt/vzp  -  tmpi*areaa*areab* 

&  resp/prwt/vzp  -  tmpi*areaa*areab*y (1) *y (1) /vzp/vzp 

c41  =  tmpi*areaa*areaa*za2/vzp/prwt 
c42  =  tmpi*areaa*areaba/vzp/prvit 
c4  =  c4]+c42 

c4  =  delta*tmpi*areaa*areaa*za2/vzp/prwt  +  delta*tmpi*areaa* 

&  areaba/vzp/prwt 

c5  =  tinpi*areaa*areaa*zaO/vzp/prwt 

11  =  zal+fk+a3*qlzp+c3*q3zp+bt*q2zp+hl*q4zp+a jl*q5zp+akl 

12  =  l-za2-a4*qlzp  -  c4*q3zp 

13  =■  zaO  +  a5*qlzp  +c5*q3zp 

bl  =  (a3*q7za+a3*q9zp+bt*q6zp+zal* (vzp-vza) +f k* (vzp-vza) 

*  +c3*q8zp+hl*qlzp+a jl*q3zp+akl* (vzp-vza) ) /vzp 
b2  =  (a4*q7za+a4*q9zp+za2* (vzp-vza) +c4*q8zp) /vzp 
b3  =  (vza+a5*q7za+a5*q9zp+za0* (vzp-vza) +c5*q8zp) /vzp 
calculate  base  pressure 

pbase= (pmean/b3  -  bl/b3  +  11/13) /(  12/13  +  b2/b3) 
calculate  breech  pressure 
pbrch=pmean/b3  -  bl/b3  -  pbase*b2/b3 
pza=zal  +  za2*pbase  +  zaO*pbrch 
write (6,7)y(3),z(l),y(l),y(2) , pmean, phase, pbrch 
write(6,7)  bl , b2, b3, 11, 12, 13 
write(6,7)  a3, a4 , a5, c3, c4 , c5 

z (1 ) = (areaa*pza+areaba*pbase-areab*resp) /prwt 
go  to  615 

use  lagrange  pressure  gradient  equation 
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if (iswl . ne . 0) go  to  600 
if (pmean . It . resp) resp  =  pmean 


calculate  base  pressure 
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tmp2  =1.0+  tmpi  /  2 . 0  /  prwt 

tnir2  =  1.0  +  tmpi  /  2.0  /  rcwt 

tmr3  =  1.0  +  tmpi  /  3.0  /  rcwt 

tmr4  =  rfor  /  areab  +  resp  -  pgas  *  air 

phase  =  pmean  /  tmr3 

1  -  tmpi  I  2.0  I  tmr2  *  (  tmr4  /  rcwt  -  resp  /  prwt) 

2  +  tmpi  /  3.0  /  tmr3  *  (  tmr4  /  rcwt  -  resp  /  prwt  /  2.0) 
phase  =  phase  /  {  tmp2  /  tmr2  -  tmpi  /  tmr3  /  prwt  /  6.0) 

c 

if (phase . gt . resp  +  l.)iswl  =  1 
c 

c  calculate  hreech  pessure 

c 

phrch  =  phase  *  tmp2  /  tmr2  +  tmpi  /  2 . 0  /  tmr2  * 

1  (tmr4  /  rcwt  -  resp  /  prwt) 

c 

c  calculate  projectile  acceleration 
c 

610  z(l)  =  areah  *  (phase  -  resp)  /  prwt 

615  if (z(l) .lt.0.0)go  to  620 

go  to  630 

620  if (iswl .eq. 0) z (1)  =  0.0 

630  if (y(l) .It.O.O)y(l)  -0.0 

z(2)  =  y(l) 
c 

c  get  burning  rate 

c 

640  do  670  m  =  1,  nprop 

z (ibrp  +  m)  -0.0 
d2xdt2(m)  =0.0 
if (ibo (m) .eq. 1)  goto  670 
do  650  )c  =  1,  nbr  (m) 

if (pmean . gt .pres (m,  k))go  to  650 
go  to  660 

650  continue 

)<  =  nbr  (m) 

660  pmix  =  pmean 

if  (igrad.eq.  3)pmix  =  phrch  -  (a)cll  *  phase  +  akl2)  /  6.  * 
1  zb  *  zb 

if (igrad.eq. 4)pmix  =  phrch  +  (alt  +  a2t  *  phase)  *  j3zb  / 
1  vzb  +  bt  *  j4zb  /  vzb 

if (pmix . It .. 99  *  pmean) pmix  =  pmean 

z(ibrp  +  m)  =  beta  (m,  )c)  *  (pmix  *  l.e  -  6)  **  alpha  (m,  k) 
abr (m)  =  alpha (m,  k) 
bbr (m)  =  beta (m,  k) 

d2xdt2(m)  =  beta (m,  k)  *  alpha (m,  k)  *  (pmix  *  l.e  -  6)  ** 
1  (alpha (m,  k)  -  1.)  *  dpmdt  *  l.e  -  6 

670  continue 

do  680  i  =  1,  nde 

d(i)  =  (z(i)  -  b(j)  +  p(i))  *  a(j) 

y(i)  =  deltat  *  d(i)  +  y(i) 

p(i)  =  3.  *  d(i)  -  ak(j)  *  z(i)  +  p(i) 

680  continue 

690  continue 
nzp=nzp+l 

if (igrad. ne . 6) goto  2003 
if (nzp. ne .nzpi) goto  2003 
dzp=zp/50 . 
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open (unit=ll, file-' dist .dat' ) 
open (unit =12, file-' press . dat' ) 
do  2001  i=l,50 
ddzp-i*dzp 

call  jint (btd, btl, zp, ddzp, nchpts, chdist,  chdiam,  bint  ,bvol) 
if (ddzp . It . za) then 

pz-pbrch-tmpi*areab*z (1 ) *bint (1) /vzp/vzp 
&  +tal*bint (1) +bt*bint (2) 

pzl-pbrch+ (a3+a5*pbrch+a4*pbase) *bint (1) +bt*bint (2) 
else 

pz-pza  +f k-tmpi*areab*z (1) *bint (10) /vzp**2 
&+tmpi*areaa*z (1) *bint (5) /vzp+bt*bint (2) +akl 
&  +hl*bint (6)-bt*2.*bint (10)-hl*bint (5)+ajl*bint (7) 
pzl=za0*pbrch+zal+za2*pbase+f k+a3*bint (10 ) 

&+a4*pbase*bint (10) +a5*pbrch*bint (10) 

&+c3*bint (5) +c4*pbase*bint (5) +c5*pbrch*bint (5) 

&+bt*bint (2) +hl*bint (6) +a jl*bint (7) +akl 
endif 

write (6, *) ' za  ' , za, '  ddzp  ',ddzp, '  pza  ',pza 
write(6,*)'  pz  ',pz,'  pzl  ',pzl 
write (11, *) ddzp 
write (12,*)pz/l.e6 
2001  continue 
2003  if (prwtO . ne .prwt) then 

if (mode .eq. 1) write (6,  1450) prwt 
i f (mode . eq . 2 ) then 

argl  =  prwt  /  0.45359237 
write (6,  1450)argl 
endif 

prwtO  -  prwt 
lines  =  lines  +  1 
endif 

t  =  t  +  deltat 
told  =  y(3) 

if(ihl.eq.l  .and.  pmean .gt .burstp) then 
write(6,  1440) 
ihl  =  2 

lines  =  lines  +  1 
endif 

if (pmaxm. gt .pmean) go  to  700 
pmaxm  =  pmean 
tpmaxm  =  y(3) 

700  if (pmaxba .gt .phase) go  to  710 
pmaxba  =  phase 
tpmxba  -  y ( 3 ) 

710  if (pmaxbr .gt .pbrch) go  to  720 
pmaxbr  -  pbrch 
tpmxbr  =  y(3) 

720  continue 

if (y (3) . It .ptime) go  to  730 
ptime  =  ptime  +  deltap 
pjt  =  y(2)  +  y(7) 
argO  -  y(3)  *  1000. 

if (mode .eq. 1) write (6,  1270)arg0,  z(l),  y(l),  pjt,  pmean,  phase, 
1  pbrch 

if (nzp. ne .nzpi) goto  2004 

write (6, *) ' dpza  ',dpza,'  pza  ',pza,'  dl  ',dl,'  gl  ',gl 
c  write (6, *) areaa*pza, areaba*pbase, areaa*pza+areaba*pbase 
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write (6, * ) ' areaa  ',areaa,'  areab  ',areab, '  areaba  ',areaba 

write(6,^)'  vza  ',vza,'  vzp  ',vzp 

write (6, *) ' qlza-qlOza' 

write  (6,  )  qlza,  q2za,  q3za,  q4za 

write  (6,  *)  qSza,  q6za,q7za,  qSza 

write (6, *) q9za, qlOza 

write  (6,*)'  qlzp-q9zp' 

write (6, * ) qlzp, q2zp, q3zp, q4zp 

write (6, *) qSzp, q6zp, q7zp, q8zp 

write (6, *)q9zp 

write  (6, *) 

write (6, *) ' zao, zal, za2, resp  ' , zaO, zal, za2,resp 
write(6, *) 'bl,b2,b3,ll  ' , bl, b2, b3, 11 
write (6,*)'12,13,a3/a4  ',12,13,a3,a4 
write(6,*)'  a5,c3,c4,c5  ' , a5, c3, c4, c5 
write (6, *) ' bt , hi, akl, a jl, fk  ' ,bt, hi, akl, a jl, fk 
2004  if (mode . eq. 2) then 

argl  =  y(l)  /  0.0254 
arg2  =•  pjt  /  0.0254 
arg3  =  pmean  /  6894.757 
arg4  =  phase  /  6894.757 
arg5  =  pbrch  /  6894.757 
arg6  =  z(l)  /  0.0254 

write  (6,  1270)arg0,  arg6,  argl,  arg2,  arg3,  arg4,  arg5 
endif 

lines  =  lines  +  1 
if (igrad. gt . 2) then 
pjt  =  y(2)  +  y(7) 
prt  =  y(10)  +  y(7) 

if (mode .eq. 1) write (6,  1280)prt,  pjt 
if (mode .eq. 2) then 

argl  =  prt  /  3.28083 
arg2  =  pjt  /  3.28083 
write(6,  1280)argl,  arg2 
endif 

lines  =  lines  +  1 
endif 

730  continue 

if (lines. gt .linmax) then 

write(6,  1236)  title, vsn 
write(6,  1230) 
if (mode .eq. 1) write (6,  1232) 
if (mode . eq. 2) write (6,  1234) 
lines  =  4 
endif 

if (t . gt . tstop) goto  740 

if(y(2)  +  y  (7)  .gt .travp)go  to  740 

rmvelo  =  y(l) 

tmvelo  =  y(3) 

rcvelo  =  y(6) 

disto  =  y(2)  +  y(7) 

go  to  250 

740  if (lines . gt . linmax-nprop-25)  write(6,  1236)  title, vsn 
write(6,  1290)t,  y(3) 

if (mode .eg. 1)  write  (6,  1300)  .maxm,  tpmaxm,  pmaxba,  tpmxba,  pmaxbr, 
1  tpmxbr 

if (mode .eq. 2) then 

pmaxm  =  pmaxm  /  6894.757 
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pmaxba  =  pmaxba  /  6894.757 
pmaxbr  =  pmaxbr  /  6894.757 

write (6,  1310 ) pmaxm,  tpmaxm,  pwnaxba,  tpmxba,  pmaxbr,  tpmxbr 
endif 

if{y(2)  +  y (7) . le.travp)goto  750 

dfract  »  (travp  -  disto)  /  (y(2)  +  y(7)  -  disto) 

rmvel  =  (y(l)  -  rmvelo)  *  dfract  +  rmvelo 

tmvel  =  (y(3)  -  tmvelo)  *  dfract  +  tmvelo 

revel  -  (y(6)  -  revelo)  *  dfract  +  revelo 

if (mode .eq. 1) write (6,  1320) rmvel,  tmvel,  revel 

i f (mode . eq . 2 ) then 

rmvel  -  rmvel  *  3.28083 
revel  =  revel  *  3.28083 
write  (6,  1330) rmvel,  tmvel,  revel 
endif 
go  to  760 

750  if (mode .eq. 1) write (6,  1340)y(l),  y(3) 
if (mode .eq. 2) then 

argl  =  y(l)  *  3.28083 
write(6,  1350)argl,  y(3) 
endif 

760  efi  =  chwi  *  forcig  /  (gamai  -  1.) 
efp  =  0.0 

do  770  i  =  1,  nprop 

efp  =  efp  +  chwp(i)  *  forcp(i)  /  (gamap(i)  -  1.0) 

770  continue 

tenerg  =  efi  +  efp 
if (mode . eq. 1) write (6,  1360) tenerg 
if (mode. eq. 2) tenerg  -  tenerg  /  0.1129848 
if (mode. eq. 2) write (6,  1370)tenerg 

tengas  ="  chwi  *  forcig  *  tgas  /  (gamai  -  1.)  /  tempi 
do  780  i  =“  1,  nprop 

tengas  =  (frac(i)  *  chwp(i)  *  forcp(i)  *  tgas  /  tempp(i)  / 

1  (gamap(i)  -  1.))  +  tengas 

780  continue 

write(6,  1380)  (i,  frac(i),  tbo(i),  i  =  1,  nprop) 
if (mode . eq. 1) write (6,  1390) 
if (mode .eq. 2) then 

tengas  =  tengas  /  0.1129848 
elpt  =  elpt  /  0.1129848 
elpr  =  elpr  /  0.1129848 
elgpm  =  elgpm  /  0.1129848 
elbr  =  elbr  /  0.1129848 
elrc  =  elrc  /  0.1129848 
elht  =  elht  /  0.1129848 
3lar  =  elar  /  0.1129848 
write (6,  1400) 
endif 

pctenl  =  tengas  /  tenerg  *  100.0 
pcten2  =  elpt  /  tenerg  *  100.0 

pcten3  =  elpr  /  tenerg  *  100.0 

pcten4  =  elgpm  /  tenerg  *  100.0 
pcten5  =  elbr  /  tenerg  *  100.0 

pcten6  =  elrc  /  tenerg  *  100.0 

pcten7  =  elht  /  tenerg  *  100.0 

pcten8  =  elar  /  tenerg  *  100.0 

write  (6,  1410) tengas,  pctenl,  elpt,  pcten2,  elpr,  pcten3, 

1  elgpm,  pcten4,  elbr,  pcten5,  elrc,  pcten6,  elht,  pcten7. 
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2  elar,  pctenS 
stop 

790  write(  *,  1420) 
stop 

800  write(  *,  1430) 

810  continue 
820  continue 
stop 

830  format  ('  input  name  of  data  file  to  be  used  as  input 
840  format  (alO) 

850  format  ('  input  name  of  output  file  ') 

860  format  (/'  the  input  file  is  ',al0/) 

870  format  ('  input  data  units  -  "metric”  or  "english"' ) 

880  format  ('  must  use  quoted  "m"  or  "e"  as  first  character  of', 

1  '  input  file'/'  to  specify  metric  or  english  input  units') 

885  format  (15a4) 

890  format  (lx, 'using  lagrange  pressure  gradient') 

900  format  (lx, 'using  chambrage  pressure  gradient') 

910  format  (lx, 'using  rga  gradient') 

920  format  (///,'  chamber  distance  m  chamber  diameter  m' , / 

1  (5x,el4 . 6, 5x,el4 . 6) ) 

925  format  (///,'  chamber  distance  in  chamber  diameter  in',/ 

1  (5x,el4 . 6, 6x, el4 . 6) ) 

930  format  (lx, 'use  first  5  points') 

940  format  (lx, '  #  points  ?  ' ) 

950  format  (lx, 'using  2  phase  gradient  equation') 

960  format  (/'  chamber  volume  m**3  ',el6.6/ 

1  '  groove  diam  m  ',el6.6/'  land  diam  m  ',el6.6/ 

1  '  groove/land  ratio  ',el6.6/'  twist  turns/caliber  ',el6.6/ 

1  '  projectile  travel  m',el6.6//'  gradient  #  ',i3,// 

1  '  variable  mass  switch', i3/'  container  switch', i7/ 

1  '  friction  factor  ',el8.6/) 

970  format  (/'  chamber  volume  in**3  ',el6.6/ 

1  '  groove  diam  in  ',el6.6/'  land  diam  in  ',el6.6/ 

1  '  groove/land  ratio  ',el6.6/'  twist  turns/caliber  ',el6.6/ 

1  '  projectile  travel  in',el6.6//'  gradient  #  ',i3,// 

1  '  variable  mass  switch', i3/'  container  switch', i7/ 

1  '  friction  factor  ',el8.6/) 

980  format  (lx, 'number  of  variable  projectile  mass  points  ',i2,/ 

1  lx,'  travel  (m)  projectile  mass  ()cg)'/ 

1  (Ix,el4.6,el4.6) ) 

985  format  (lx, 'number  of  variable  projectile  mass  points  ',i2,/ 

1  lx,'  travel  (in)  projectile  mass  (lb)'/ 

1  (Ix,el4.6,el4.6) ) 

990  format  (lx, ' canister  burst  pressure  (mpa) ' , el4 . 6/ 

1  lx, 'canister  volume  (m''3)  ',el4.6/ 

1  lx, 'canister  diameter  (m)  ',el4.6) 

1000  format  (lx, ' canister  burst  pressure  (psi) ' , el4 . 6/ 

1  lx, 'canister  volume  (in'^l)  ',el4.6/ 

1  lx, 'canister  diameter  (in)  ',el4.6) 

1010  format  (/'  projectile  mass  )cg' ,  34x,  el 4 . 6/ 

1  '  switch  to  calculate  energy  lost  to  air  resistance  ',i3/ 

1  '  fraction  of  work  against  bore  used  to  heat  the  tube',el4.6/ 

1  '  gas  pressure  mpa  ',el4.6) 

1020  format  (/'  projectile  mass  lb' , 34x, el4 . 6/ 

1  '  switch  to  calculate  energy  lost  to  air  resistance  ',i3/ 

1  '  fraction  of  work  against  bore  used  to  heat  the  tube',el4.6/ 

1  '  gas  pressure  psi  ',el4.6) 
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format  (/'  number  barrel  resistance  points  ' » i2 

L  '  bore  resistance  mpa  -  travel  m  ' / (3x, el4 . 6, 8 
format  (/'  number  barrel  resistance  points  ' » i2 

L  '  bore  resistance  psi  -  travel  inches  '/(3x,el4.6 
format  (lx) 
format  {/ 

L  '  mass  of  recoiling  parts  )cg  ',el4.6/ 

L  '  number  of  recoil  point  pairs  ',i3/ 

L  '  recoil  force  J' , '  recoil  time  sec' /, (lx, el4 . 6, 3: 
format  (/ 

L  '  mass  of  recoiling  parts  lb  ',el4.6/ 

L  '  number  of  recoil  point  pairs  ',i3/ 

L  '  recoil  force  in-lb'/'  recoil  time  sec'/ 

L  (lx, el4 . 6, 3x,el4 .6) ) 
format  (/ 

L  '  free  convective  heat  transfer  coefficient  w/m^2  k 
L  '  chamber  wall  thickness  m 

L  '  heat  capacity  of  steel  of  cheimber  wall  j/kg  k 
L  '  initial  temperature  of  chamber  wall  k 
L  '  heat  loss  coefficient 
L  '  density  of  chamber  wall  steel  kg/m'"! 
format  (/ 

L  '  free  convective  heat  transfer  coef  in-lb/in''2-s-k 
L  '  chamber  wall  thickness  (inches) 

L  '  heat  capacity  of  steel  of  chamber  wall  in-lb/lb-k 
L  '  initial  temperature  of  chamber  wall  k 

L  '  heat  loss  coefficient 

L  '  density  of  chamber  wall  steel  lb/in^3 

format  ( 

L  '  impetus  of  igniter  propellant  j/kg 

L  '  adiabatic  flame  temperature  of  igniter  propellant 
L  '  covolume  of  igniter  m'"3/kg 

L  '  ratio  of  specific  heats  for  igniter 

L  '  initial  mass  of  igniter  kg 

format  ( 

L  '  impetus  of  igniter  propellant  ft-lb/lb 

L  '  adiabatic  flame  temperature  of  igniter  propellant 
L  '  covolume  of  igniter  ft''3/lb 

L  '  ratio  of  specific  heats  for  igniter 

L  '  initial  mass  of  igniter  lb 

format  {/'  there  are  ',i2,'  propellants'//) 
format  (('  for  propellant  number', i2// 

L  '  impetus  of  propellant  j/kg  ',el4.6/ 

L  '  adiabatic  temperature  of  propellant  k  ',el4.6/ 

L  '  covolume  of  propellant  m''3/kg  ',el4.6/ 

L  '  ratio  of  specific  heats  for  propellant  ',el4.6/ 

L  '  initial  mass  of  propellant  kg  ',el4.6/ 

L  '  density  of  propellant  kg/m''3  ',el4.6/ 

L  '  number  of  perforations  of  propellant  ',13/ 

L  '  length  of  propellant  grain  m  ',el4.6/ 

L  '  perforation  diameter  m  ',el4.6/ 

L  '  outside  diameter  of  propellant  grain  m  ',el4.6/) 

format  (('  for  propellant  number',i2// 

L  '  impetus  of  propellant  ft-lb/lb  ',el4.6/ 

L  '  adiabatic  temperature  of  propellant  k  ',el4.6/ 

L  '  covolume  of  propellant  in''3/lb  ',el4.6/ 

L  '  ratio  of  specific  heats  for  propellant  ',el4.6/ 

L  '  initial  mass  of  propellant  lb  ',el4.6/ 


X, el4 . 6 ) 
/ 

,e22.6) ] 


X, el4 . 6) 


'  ,el4 
'  ,el4 
'  ,el4 
'  ,el4 
'  ,el4 
' ,el4 . 

' ,el4 , 
' ,el4 , 
'  ,el4 
'  ,el4 
'  ,el4 
',el4, 

'  ,el4 
k' ,el4 
'  ,el4 
'  ,el4 
',el4, 

' ,el4 , 
k' ,el4 
'  ,el4 
'  ,el4 
',el4. 
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1  '  density  of  propellant  lb/in''3  ',el4.6/ 

1  '  number  of  perforations  of  propellant  ',13/ 

1  '  length  of  propellant  grain  in  ',el4.6/ 

1  '  perforation  diameter  in  ',el4.6/ 

1  '  outside  diameter  of  propellant  grain  in',el4.6/)M 

1150  format  (/'  for  propellant  ',i2,'  the  total  number  of  ^  ains' 

1  is  ',el4.6) 

1160  format  ('  number  of  burning  rate  points', 12/ 

1  '  exponent  coefficient  pressure' / 

1  '  -  m/sec-mpa''ai  mpa' ) 

1170  format  ('  number  of  burning  rate  points', 12/ 

1  '  exponent  coefficient  pressure'/ 

1  '  -  in/sec-psi''ai  psi') 

1180  format  (lx, el4 . 6, 5x, el4 . 6, 4x, el4 . 6) 

1190  format  {'  time  increment  msec  ',el4.6/ 

1  '  print  increment  msec  ',el4.6/ 

1  '  time  to  stop  calculation  msec',el4.6) 

1200  format  (lx, 'end  input  data  —  i.b.  calculation  start') 

1210  format  (/'  area  bore  m^2  ',e27.6/'  pressure  from  ign  pa',e21.6/ 
1  ,'  volume  of  unburnt  prop  m^3  ',el4.6/ 

1  '  init  cham  vol-cov  ign  m''3  ',el5.6) 

1220  format  (/'  area  bore  in''2  ',e29.6/'  pressure  from  ign  psi',e23.6/ 
1  ,'  volume  of  unburnt  prop  in''3  ',el6.6/ 

1  '  init  cham  vol-cov  ign  in^3  ',el7.6) 

1230  format  {/'  time  accel  velocity  distance  pr  mean', 

1  '  pr  base  pr  brch' ) 

1232  format  {  '  (ms)  (m/s/s)  (m/s)  (m)  (Pa)  ', 

1  '  (Pa)  (Pa)') 

1234  format  (  '  (ms)  (in/s/s)  (in/s)  (in)  (psi)  ', 

1  '  (psi)  (psi)') 

1236  format  ( Ihl , 3x, 15a4 , '  rga.',a4) 

1240  format  ('  propellant' , i2, '  has  slivered') 

1250  format  ('  propellant' , i2, '  has  burned  out') 

1270  format  ( lx, lp7el 1 . 4 ) 

1280  format  (lx, 'prop  travel ', ell . 4 ,' pro j  travel ' , el  1 . 4 ) 

1290  format  (/lx,'  deltat  t',  el4.6,  '  intg  t',el4.6) 

1300  format  (/'  pmaxmean  pa  ',lpel4.7,'  time  at  pmaxmean  sec  ', 

1  0pel6.6/'  pmaxbase  pa  ',lpel4.7,'  time  at  pmaxbase  sec  ', 

1  Opel 6. 6/'  pmaxbreech  pa  ',lpel4.7,'  time  at  pmaxbreech  sec  ', 

1  Opel  4.6) 

1310  format  (/'  pmaxmean  psi', f 14. 3,'  time  at  pmaxmean  sec  ', 

1  0pel6.6/'  pmaxbase  psi', f 14. 3,'  time  at  pmaxbase  sec  ', 

1  0pel6.6/'  pmaxbreech  psi', f 14. 3,'  time  at  pmaxbreech  sec  ', 

1  0pel4 . 6) 

1320  format  (/lx, 'muzzle  velocity  m/s  ',el4.6,'  time  of  muzzle  exit  ', 
1  el4.6,'  sec' // lx, ' recoil  velocity  m/s  ',el4.6) 

1330  format  (/lx, 'muzzle  velocity  ft/s  ',el4.6,'  time  of  muzzle  exit  ', 
1  el4.6,'  sec' // lx, ' recoil  velocity  ft/s  ',el4.6) 

1340  format  (/  '  velocity  of  projectile  m/s  ',el4.6,'  at  time  ',el4.6, 
1  '  sec' ) 

1350  format  (/  '  velocity  of  projectile  ft/s  ',el4.6,'  at  time  ',el4.6, 
1  '  sec' ) 

1360  format  (/lx, 'total  initial  energy  available  j  =  ',el4.6/) 

1370  format  (/lx, 'total  initial  energy  available  in-lb  =  ',el4.6/) 

1380  format  ('  propellant  mass  fraction  burnt  time  <sec)'/ 

1  (4x,  i2,  9x,  el4 . 6,  5x,  el4 . 6)  ) 

1390  format  (/'  **  energy  summary  **',  23x,  ' joules' ,  1 lx,  '  °  '  ) 

1400  format  (/'  **  energy  summary  **',23x,'  in-lb' , 1  lx ,'%' ) 
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1410  format  (lx, 'total  energy  remaining  in  gas',llx,' 

1  /lx, 'energy  loss  from  pro3eci-ile  translation 

1  /lx, 'energy  loss  frcr  projectile  rotation 

1  /lx, 'energy  lost  to  gas  and  propellant  motion 
1  /lx, 'energy  lost  to  bore  resistance 
1  /lx, 'energy  lost  to  recoil 
1  /!>., 'energy  loss  from  heat  transfer 

1  /lx, 'energy  lost  to  air  resistance 
1420  format  (]x,'end  of  file  encounter') 

1430  format  (lx, 'read  or  write  error') 

1440  format  ('  canister  burst  pressure  achieved') 

1450  forr at  ('  projectile  mass  transition  point  -  new  mass  =  ', 

1  lpell.4) 
end 

subroutine  prf710(pd,  gd,  ql,  np,  x,  frac,  vol,  surf,  dsdx) 
common  nsl,  kpr,  fracsl(lO),  dsdxsl(lO),  surfsi(lO),  nslpIlO), 

1  tsl(lO),  pbrch,  phase,  pmean,  bbr(lO),  abr(lO),  deltat,  yar(20), 

1  igrad 

dimension  ts(lO),  coef(lO) 
pi  =  3.141593 
nsl  =  0 
c 

c  pd=perf orat ion  diameter 

c  gd=outer  dia 

c  gl=grain  length 

c  np=number  of  perfs 

c 

c  surf=output  surfacf;  area 

c  frac*output  mass  fraction  of  propellant  burned 

c 

c  w  =  web  =  distance  between  perforation  edges 

c  =  distance  between  outside  perf  edge  and  edge  of  grain 

c 

c  p  =  distance  between  perforation  centers 

c 

c  xl  =  distance  to  inner  sliver  burnout 
c 

c  x2  =  distance  to  outer  sliver  burnout  (frac=l) 

c 

if(np.eq.O)  go  to  70 
if (np.eq. 1) go  to  90 
if (np . eq. 2) go  to  210 
if (np.eq. 7) go  to  10 
if (np.eq. 19) go  to  110 
if (np.eq. 15)go  to  180 
write(6,  220) 
stop 

10  w  =  (gd  -  3.  *  pd)  /  4. 
d  =  w  +  pd 
sqr3  =  sqrt  (3 . ) 
xl  =  d  /  sqr3  -  pd  /  2 . 

x2  =  (14.  -  3.  *  sqr3)  *  d  /  13.  -  pd  /  2. 
vO  =  pi  /  4 .  *  gl  *  (gd  *  gd  -  7.  *  pd  *  pd) 
sO  =  2.  *  vO  /  gl  +  pi  *  gl  *  (gd  +7.  *  pd) 
if  (x.gt.w  /  2.  +  .0000001)  goto  20 

vol  =  pi  /  4 .  *  (gl  -  2.  *  xj  *  ( (gd  -2.  *x)  **2-7.  *  (pd  + 
1  2.  *  X)  **  2) 

surf  =2.  *  vol  /  (gl  -  2.  *  x)  +  pi  *  (gl  -  2.  *  x)  *  ( (gd  -  2. 


=  ',el4.6,fll.4 
=  ' , el4 . 6,  fll .  4 
»  ',el4.6,fll.4 
=  ',el4.6,fll.4 
-  ',el4.6,fll.4 
=  ',el4.6,fll.4 
“  ',el4.6,fll.4 
=  ' ,el4  .6, fll. 4) 
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gl  +  18.  *  X) 


I 

i 

I 

I 


1  *  X)  +7.  *  (pd  +  2.  *  X)) 
frac  »  1.  -  vol  /  vO 

dsdx  --4*pi*  (gd+7.  *pd-3.  * 
dsdxsl(kpr)  =  dsdx 
fracsl(kpr)  =  frac 
surfsl(kpr)  =  surf 
return 
20  nsl  =  1 

coef(kpr)  =  0. 

if (igrad.eq. 1 .or . igrad.eq.2)go  to  40 
if (nslp (kpr ) . eq. 1 ) goto  30 
tsl (kpr)  =  yar (3) 
ts(kpr)  =w/2.  *(-l.  +  (pbrch  /  pmean)  **  abr(kpr))  / 

1  (bbr(kpr)  *  (phase  *  l.e  -  6)  **  abr(kpr)) 

30  continue 

coef(kpr)  =  (ts(kpr)  +  tsl (kpr)  -  (deltat  +  yar (3)))  /  ts(kpr) 
if (coef (kpr) .gt . 1 . ) coef (kpr)  =  1. 
if (coef (kpr ) . It . 0 . ) coef (kpr)  =  0. 

40  if (x . ge . x2 ) goto  60 
si  =  0. 
s2  =  0. 
vl  =  0. 
v2  =  0. 
dsldx  =  0. 
ds2dx  =  0 . 

y  =  sqrt ( (pd  +2.  *  x)  **  2  -  d  *  d) 
theta  =  atan(y  /  d) 


al  = 

theta  /  4 . 

*  (pd  +  2 . 

*  x)  ** 

2  -  d 

/  4. 

*  y 

if  (X, 

.  ge . xl) goto 

50 

vl  = 

3.  /  4.  * 

(gl  - 

2.  * 

X) 

vl  = 

vl  *  (2.  * 

sqr3 

*  d  * 

d  -  pi  * 

(pd  + 

2.  * 

X) 

★  ii 

'  2 

+  24. 

*  al) 

si  = 

2.  *  vl  / 

(gl  - 

2.  * 

X) 

si  = 

si  +  3.  * 

(gl  - 

2.  * 

x)  *  (pi 

-  6.  * 

theta) 

★ 

(pd 

+  2. 

*  X) 

yl  = 

sqrt ( (gd  - 

2.  + 

x)  ** 

2  -  (5.  ’ 

'  d  -  2 

★ 

(pd 

+ 

2. 

*  X)) 

**  2) 

chi  = 

=  atan(yl  / 

(5. 

*  d  - 

2.  *  (pd 

+  2.  * 

X)  ) 

) 

y2  = 

sqrt ( (pd  + 

2.  * 

x)  ** 

2  -  (3.  ^ 

►  d  -  2 

★ 

(pd 

+ 

2. 

*  X)) 

**  2) 

phi  = 

=  atan(y2  / 

(3. 

*  d  - 

2 .  *  (pd 

+  2.  * 

X)) 

) 

a2  = 

phi  *  (pd 

+  2. 

*  X)  * 

*  2  -  chi 

*  (gd 

-  2 

★ 

X) 

★  ★ 

2 

a2  = 

(a2  +  2.  * 

sqr3 

*  d  * 

sqrt ( ( 3 . 

*  d  - 

pd 

-  2. 

★ 

X) 

*  (3. 

,  *  d 

1  -  gd  +  2 .  *  x) 

))  / 

8. 

v2  = 

pi  *  (gd  - 

2.  * 

X)  ** 

2-6.  * 

sqr3 

*  d 

*  d 

- 

4  . 

*  pi  ^ 

'  (pd 

1  +  2. 

,  *  X)  **  2 

v2  = 

(v2  +  24. 

*  (al 

+  2. 

*  a2) )  * 

(gl  - 

2.  * 

X) 

/ 

4  . 

s2  = 

2.  *  v2  / 

(gl  - 

2.  * 

X) 

s2  = 

s2  +  (gl  - 

2.  * 

X)  * 

( (pi  -  6 . 

*  chi 

)  * 

(gd 

- 

2  . 

*  X)  -1 

-  2.  * 

1  (2. 

*  pi  -  3 . 

*  phi 

-  3. 

*  theta  ) 

*  (pd 

+  2 

★ 

X) 

) 

vol  =  vl  +  v2 

surf  =  si  +  s2 

frac  =  1.  -  vol  /  vO 

dsdx  =  -  surf  /  (x2  -  x) 


dsdx  =  coef (kpr)  * 
dsdxsl(kpr)  =  dsdx 

dsdxsl(kpr)  + 

(1. 

-  coef (kpr) ) 

*  dsdx 

frac  =  coef (kpr)  * 
fracsl(kpr)  =  frac 

fracsl (kpr)  + 

(1. 

-  coef (kpr) ) 

*  frac 

surf  =  coef (kpr)  * 

surfsl(kpr)  + 

(1. 

-  coef (kpr) ) 

*  surf 
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surf 


surfsl(kpr)  = 
return 
60  vol  =  0 . 
surf  =  0. 

frac  =  fracsl(kpr)  *  coef(kpr)  +  1.  -  coef{kpr) 

fracsl{kpr)  =  frac 

if (frac .gt .. 9999)  frac  =  1. 

if (frac. gt. .9999)return 

dsdx  =  0 . 

dsdx  =  dsdxsl(kpr)  *  coef (kpr) 
dsdxsl(kpr)  =  dsdx 
if (abs (dsdx) . It . 1 . ) dsdx  =  0. 
surf  =  surf si (kpr)  *  coef (kpr) 
surfsl(kpr)  =  surf 
return 
c 

c  zero  perf  calculations  start  here, 
c 

70  if (gd  -  2.  *  x.le.0.0)  go  to  80 
vO  =  pi  *  gd  *  gd  /  4 .  *  gl 

vol  =  pi  *  (gd  -  2.  *  x)  **  2  /  4.  *  (gl  -  2.  *  x) 
frac  =  1 .  -  vol  /  vO 

surf  =  pi  /  2 .  *  (gd  -  2.  *  x)  **  2  +  pi  *  (gd  -  2.  *  x)  *  (gl  - 
1  2.  *  X) 

dsdx  =  -  2.  *  pi  *  (gd  +  gl  -  6 .  *  x) 
return 
c 

c  grain  completely  burned 

c 

80  surf  «  0. 
frac  =  1.0 
vol  =  0 . 
dsdx  =  0 . 
nsl  =  1 
return 
c 

c  one  perf  calculation  starts  here 

c 

90  if (gd  -  pd  -  4.  *  x.le.0.0)  goto  80 

vO  =  pi  /  4 .  *  (gd  *  gd  -  pd  *  pd)  *  gl 

vol  =  pi  /  4.  *  ((gd  -  2.  *  x)  **  2  -  (pd  +  2 .  *  x)  **  2)  *  (gl 
1  2.  *  X) 

frac  =  1.  -  vol  /  vO 

surf  =  pi  /  2 .  *  ( (gd  -  2.  *  x)  **  2  -  (pd  +2.  *  x)  **  2) 

surf  =  surf  +  pi  *  (gd  -  2.  *  x)  *  (gl  -  2.  *  x) 

surf  =  surf  +  pi  *  (pd  +2.  *  x)  *  (gl  -  2.  *  x) 

dsdx  =  -  4 .  *  pi  *  (gd  +  pd) 

return 
c 

c  below  is  the  calculation  for  the  cylindrical  19  perf  grain, 
c  programmer : f red  robbins 

c  input 

c 

c  p  =  perf  diameter 

c  d  =  grain  diameter 

c  gl  “  grain  length 

c  X  =  distance  burnt 

c 


90 


c 

c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


c 

c 

c 

c 

c 

c 

c 


c 

c 

c 


output 

vol  “  the  volume  of  one  grain  at  x. 

surf  =  the  surface  area  of  one  grain  at  x. 

frac  =  the  fraction  of  grain  burnt  at  x. 

w=web 

110  p  =  pd 
d  =  gd 

w=  (d-5.  *p)  /  6. 
pi  =  3.141592654 
sqrt3  =  sqrt (3 . ) 
sqrt5  =  sqrt (5 . ) 
sqrt6  =  sqrt ( 6 . ) 
sqrtlO  =  sqrt (10.) 

initial  volume  and  surface  area 

vO  =  pi  /  4 .  *  gl  *  (d  *  d  -  19.  *  p  *  p) 

sO  =  2.  *  vO  /  gl  +  pi  *  gl  *  (d  +  19.  *  p) 

xl  =  distance  to  inner  sliverr  burnout 

x2  =  distance  to  outer  sliver  burnout 

dbc  =  distance  between  perforation  centers 

assumes  burnout  does  not  occur  in  longitudinal  direction 

wl  =  secondary  web 

dbc  =  w  +  p 

wl=0.5*  (d-p-2.  *  sqrt3  *  dbc) 
xl  =  dbc  /  sqrt 3  -  p  /  2. 

x2  =  0.25  *  (dbc  *  (6.  -  sqrtlO)  -  2.  *  p) 

if(x.gt.w  /  2 . ) go  to  120 

not  slivered  yet 


vol  = 

=  pi  / 

4.  * 

(gl  '  2. 

★ 

X) 

★ 

((d  - 

2.  * 

X)  **  2  -  19 

1  *  X) 

**  2) 

surf 

=  2.  * 

vol 

/ 

(gl  - 

2. 

★ 

X) 

+  pi 

*  (gl 

-  2.  *  X)  * 

1  X  + 

19.  * 

(P  + 

2  , 

.  *  X)  1 

1 

dsdx 

=  pi  * 

{  - 

4 

*  d  + 

36 

* 

gl 

-  76 

*  p  - 

216  *  X) 

frac 

=  1.  - 

vol 

/ 

vO 

dsdxsl(kpr)  =  dsdx 
fracsl(kpr)  =  frac 
surfsl(kpr)  =  surf 
return 


*  (p  +  2. 

(d  -  2.  * 


c 

c  vl=total  volume  of  inner  sliver,  v2=total  volume  of  outer  sliver 

c  sl=total  surface  area  of  inner  slivers,  s2=total  surface  area  of 

c  outer  slivers 

c 


120  vl 

=  0 

v2 

=  0 

si 

=  0 

s2 

=  0 

delta  =  0  . 
chi  =  0 . 
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abr (kpr) )  / 


nsl  =  1 

coef(kpr)  =  0. 

if (igrad.eq. 1 .or . igrad.eq.2)go  to  140 
if (nslp (kpr ) . eg. 1 ) goto  130 
tsl(kpr)  =  yar(3) 

ts(kpr)  =w/2.  *(-l.  +  (pbrch  /  pmean) 

1  (bbr(kpr)  *  (phase  *  l.e  -  6)  **  abr (kpr)) 

130  continue 

coef(kpr)  =  (ts(kpr)  +  tsl(kpr)  -  (deltat  +  yar(3)))  /  ts(kpr) 
if (coef (kpr) .gt . 1 . ) coef (kpr)  =  1. 
if (coef (kpr ) . It . 0 . ) coef (kpr)  =  0. 

140  a3  =  0. 

if (x . ge . x2) go  to  170 

theta  =  acos(dbc  /  (p  +  2 .  *  x) ) 

al  =  theta  /  4.  *  (p  +  2 .  *  x)  **  2  -  dbc  /  4.  *  sqrt ( (p  +2.  *  x) 
1  **  2  -  dbc  *  dbc) 
if (x . gt . xl) go  to  150 

vl  =  3.  *  (gl  -  2.  *  x)  *  (2.  *  sqrt3  *  dbc  *  dbc  -  pi  *  (p  +  2. 

1  *  x)  2  +  24  *  al) 

si  =  2.  *  vl  /  (gl  -  2.  *  X)  +  12.  *  (gl  -  2.  *  x)  *  (pi  -  6.  * 

1  theta)  *  (p  +  2.  *  x) 

150  phi  =  acos((5.  *  d  -  13 .  *  p  -  36 .  *  x)  /  (12.  *  (p  +  2 .  *  x) ) ) 
xi  =  acos((13.  *  d  -  5.  *  p  -  36.  *  x)  /  (12,  *  (d  -  2.  *  x) ) ) 
if(x.le.wl  /  2.)go  to  160 

delta  =  acos((2.  *d-p-6.  *x)  /  (sqrt3  *  (d  -  2.  *  x) ) ) 
chi  =  acos((d  -2.  *p-6.  *x)  /  (sqrt3  *  (p  +  2 .  *  x) ) ) 
a3  =  .125  *  (chi  *  (p  +  2 ,  *  x)  **  2  -  delta  *  (d  -  2 .  *  x)  **  2 

1  +  2.  *  sqrt6  *  dbc  *  sqrt (6.  *  dbc  *  (p  +  2.  *  x  -  dbc)  -  (p  +  2 . 

1  *  X)  **  2) ) 

160  a2  =  ,125  *  (phi  *  (p  +  2.  *  x)  **  2  -  xi  *  (d  -  2.  *  x)  **2+2. 
1  *  sqrt 5  *  dbc  *  sqrt ((5.  *  dbc  -p-2.  *x)  *  (5.  *  dbc  -  d  +  2. 
1  *  X)  )  ) 

v2  =  .25  *  (gl  -  2.  *  X)  *  (pi  *  (d  -  2 .  *  x)  *  *  2  -  7 .  *  pi  * 

1  (p+2.  *  x)  **2-24.  *  sqrt3  *  dbc  *  dbc  +48.  *  (al  +  a2  +  a3) ) 
s2  =  2.  *  v2  /  (gl  -  2.  *  x)  +  (gl  -  2 .  *  x)  *  ((d  -  2.  *  x)  *  (pi 

1  -  6.  *  (xi  +  delta))  +  (p  +  2.  *  x)  *  (7,  *  pi  -  6 .  *  (2.  *  theta 

1  +  chi  +  phi) ) ) 
vol  =  vl  +  v2 
surf  =  si  +  s2 
dsdx  =  -  surf  /  (x2  -  x) 
frac  =  1.  -  vol  /  vO 


dsdx  =  coef (kpr)  * 
dsdxsl(kpr)  =  dsdx 

dsdxsl (kpr) 

+ 

(1. 

-  coef (kpr) ) 

*  dsdx 

frac  =  coef (kpr)  * 
fracsl(kpr)  =  frac 

f racsl (kpr) 

+ 

(1. 

-  coef (kpr) ) 

*  frac 

surf  =  coef (kpr)  * 
surfsl(kpr)  =  surf 

surf si (kpr) 

+ 

(1. 

-  coef (kpr) ) 

*  surf 

return 
170  vol  =  0. 
surf  =  0. 

frac  =  fracsl(kpr)  *  coef (kpr)  +  1.  -  coef (kpr) 

fracsl(kpr)  =  frac 

if (frac .gt ,. 9999)  frac  =  1. 

if (frac. gt. .9999)return 

dsdx  =  0 . 
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dsdx  *  dsdxsl{kpr)  *  coef(kpr) 
dsdxsl(kpr)  =  dsdx 
if (abs (dsdx) . It . 1 . ) dsdx  »  0. 
surf  surfsl(kpr)  *  coef  (kpr) 
surfsl(kpr)  »  surf 
return 
c 

c  below  is  the  calculation  for  the  19  perf  hex  grain, 
c  programmer : karen  a.  cieslewicz<std. cont .> 
c 

c  translation  of  the  input  values, 
c  p=  perf  diameter 

c  d=  grain  diameter 

c  gl=  grain  length 

c  x=  distance  burnt 

c 

translation  of  the  output  values. 
vol=  volume  of  one  grain  at  x. 
c  surf=  surface  area  of  one  grain  at  x. 

c  frac=  mass  fraction  of  the  grain  burnt  at  x. 

c 

c  assignment  statement  for  pi. 

180  pi  =  3.141592654 
sqrt3  =  sqrt (3 . ) 
p  =  pd 
d  =  gd 
c 

c  d=6w  +  5p  is  the  statement  for  the  grain  diameter  which  will  be 
c  used  to  calculate  the  web. 
c 

c  to  calculate  the  web. 

w«  (d-5.  *p)  /6. 
c 

c  below  is  the  equation  to  calculate  the  distance  between  the  perf  cen- 
c  ters. 

dpc  =  p  +  w 

c  to  calculate  the  grain  diameter  between  the  flats. 

f  =  2.  *  (sqrt3  *  dpc  +  p  /  2.  +  w) 
c 

c  to  calculate  the  distance  burnt, 
xl  =  dpc  /  sqrt3  -  p  /  2. 
x2  =  0.125  *  (5.  *  dpc  -  4.  *  p) 
c 

c  to  calculate  the  area. 

a  =  sqrt3/3.  *  ( (w  +  p  /  2.)  **  2)  -  pi  /  6 .  *  ( (w  +  p  /  2.)  **2) 
c 

c  to  calculate  the  initial  volume  of  the  sharp  corner  grain. 

vs  =  gl  /  4 .  *  (2.  *  sqrt3  *  f  **  2  -  19.  *  pi  *  p  **  2) 
c 

c  to  calculate  the  volume  that  will  be  removed  from  the  grain. 

vr  =  6.  *  a  *  gl 
c 

c  to  calculate  the  initial  volume  for  the  19hex  grain  with  rounded 
c  corners. 

VO  =  vs  -  vr 
c 

c  to  calculate  the  initial  surface  area  of  the  sharp  corner  grain, 
ss  =  2.  *  vs  /  gl  +  gl  *  (2.  *  sqrt3  *  f  +  19.  *  pi  *  p) 
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c 

c  to  calculate  the  surface  area  that  will  be  removed  from  the  grain. 

sr  =  12.  *  a  +  gl  *  (w  +  p  /  2.)  *  (4.  *  sqrt3  -  2.  *  pi) 
c 

c  to  calculate  the  initial  surface  area  for  the  19hex  grian  with  rounded 
c  corners. 

so  -  ss  -  sr 
c 

c  to  calculate  the  unknows  of  the  grain  under  the  condition  x.le.5*w. 
if (O.le.x.and.x.le.w  /  2.)  then 

a  “  sqrt3  /3.  *  (w-2.  *x+  (p+2.  *x)  /  2.)  **2-pi/ 

1  6.  (w  -  2.  *  X  +  (p  +  2.  *  X)  /  2.)  **  2 

c  to  calculate  the  volume  that  will  be  removed  from  the  sharp  corner  gra 
vr  =  6.  *  a  *  (gl  -  2.  *  x) 

c  to  calculate  the  volume  for  the  sharp  corner  grain  at  some  distance  bu 
vn  =  .25  *  (gl  -  2.  *  x)  *  (2.  *  sqrt3  *  (f  -  2.  *  x)  **  2.  - 

1  19.  pi  *  (p  +  2.  *  X)  **  2.) 

c 

c  to  calculate  the  volume  for  the  19hex  grain  with  rounded  corners. 

V  =  vn  -  vr 
c 

c  to  calculate  the  surface  area  that  will  be  removed  from  the  sharp 
c  corner  grain. 

sr  =  12.  *  a  +  (gl  -2.  *x)  *  (w-2*x+  (p+2.  *x)  12.) 
1  *  (4.  *  sqrt3  -  2.  *  pi) 

c  to  calculate  the  surface  area  for  the  sharp  corner  grain. 

sn  =  2.  *  V  /  (gl  -  2.  *  x)  +  (gl  -  2.  *  x)  *  (2.  *  sqrt3  *  (f 
1  -  2.  *  x)  +  19.  *  pi  *  (p  +  2.  *  X)) 

to  calculate  the  surface  area  for  19hex  grain  with  rounded  corners, 
s  =  sn  -  sr 
c  to  calculate  the  mass  fraction, 
frac  =  1  -  V  /  VO 

dsdx  =  -  8.  *  sqrt3  *  (f  -  2.  *  x)  -  76.  *  pi  *  (p+2.  *  x) 

1  +  (gl  -2.  *x)  *  (-4.  *  sqrt3  +38.  *  pi)  +  16  *  sqrt3  * 

1  (w  +  p  /  2 .  -  x)  -  8 .  *  pi  *  (w  +  p  /  2 .  -  x)  +  (gl  -  2.  *  x) 

1  *  (4.  *  sqrt3  -  2.  *  pi) 

surf  =  s 
vol  =  V 

dsdxsl(kpr)  =  dsdx 
fracsl(kpr)  «  frac 
surfsl(kpr)  =  surf 
return 
endif 
c 

c  due  to  the  cross  section  at  the  sliver  point  x“.5*w  there  will  be  24 
c  identical  inner  slivers, 12  identical  side  slivers,  after  slivering  th 
c  surface  area  and  the  volume  function  become  more  complex.  each  type 
o 

c  sliver  will  be  treated  separately  and  later  the  volumes  will  be  combin 
c  to  complete  the  function, 
c 

c  to  calculate  the  12  identical  side  slivers  for  the  grain  x-.5/w. 
nsl  -  1 
coef(kpr)  =  0. 

if (igrad.eq. 1 .or . igrad.eq. 2) go  to  200 
if (nslp(kpr) .eq.l)goto  190 
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•k  ★ 


abr(kpr))  / 


tsl  (kpr)  =  yar (3) 

ts(kpr)  =w/2.  *  (-1.  +  (pbrch  /  pmean) 

1  (bbr(kpr)  *  (phase  *  l.e  -  6)  **  abr(kpr)) 

190  continue 

coef{kpr)  =  (ts(kpr)  +  tsl (kpr)  -  (deltat  +  yar  (3)))  /  ts(kpr) 
if (coef (kpr) . gt . 1 . ) coef (kpr)  =  1. 
if (coef (kpr) . It . 0 . ) coef (kpr)  =  0. 

200  if(w  /  2 . It . X . and. X • It . xl .and. X . It . x2)  then 
c 

c  to  calculate  the  areas  of  the  grain. 

a  =  sqrt3  /  3.  *  (w  -  2 .  *  x  +  (p  +  2 .  *  x)  /  2.)  **  2  -  pi  / 
1  6.  *  (w  -  2.  *  X  +  (p  +  2.  *  X)  /  2.)  **  2 

theta  =  acos(dpc  /  (p  +  2 .  *  x) ) 

al  =  theta  /  4.  *  (p  +  2.  *  x)  **  2  -  dpc  /  4.  *  sqrt ( (p  +  2. 
1  *  X)  **  2  -  dpc  **  2) 

omega  =  acos(2.  *  dpc  /  (p  +  2.  *  x)  -  1.) 

a2  =  0.125  *  (p  +  2.  *  x)  *  ((p  +  2.  *  x)  *  (omega  +  sin (omega 
1  ))  -  2.  *  dpc  *  sin (omega)) 

c  to  calculate  the  volumes  of  the  grain. 

vl  =  3 .  *  (gl  -  2.  *  x)  *  (2.  *  sqrt3  *  dpc  **  2  -  pi  *  (p  + 

1  2.  *  X)  **  2  +  24.  *  al) 

v2  =  6 .  *  (gl  -  2.  *  x)  *  (2.  *  dpc  **  2  -  dpc  *  (p  +  2.  *  x) 
1  -  pi  /  4.  *  (p  +  2.  *  X)  **  2  +  2.  *  al  +  4.  *  a2) 

c  to  calculate  the  surface  areas  of  the  grain. 

si  =  2.  *  vl  /  (gl  -  2.  *  x)  +  12.  *  (gl  -  2.  *  x)  *  (pi  -  6. 
1  *  theta)  *  (p  +  2.  *  x) 

s2  =  2 .  *  v2  /  (gl  -  2.  *  x)  +  12.  *  (gl  -  2.  *  x)  *  (dpc  +  (p 
1  +  2.  *  x)  *  (pi  /  2.  -  omega  -  theta  -  sin (omega))) 

c  to  calculate  the  total  volume  and  total  surface  area, 
vf  =  vl  +  v2 
sf  =  si  +  s2 

c  to  calculate  the  mass  fraction, 
frac  =  1 .  -  vf  /  VO 
surf  =  sf 

dsdx  =  -  surf  /  (x2  -  x) 
vol  =  vf 

dsdx  =  coef (kpr)  *  dsdxsl(kpr)  +  (1.  -  coef (kpr))  *  dsdx 
dsdxsl(kpr)  =  dsdx 

frac  =  coef  (kpr)  *  fracsl(kpr)  +  (1.  -  coef (kpr))  *  frac 
fracsl(kpr)  =  frac 

surf  =  coef (kpr)  *  surf si (kpr)  +  (1.  -  coef (kpr))  *  surf 
surf si (kpr)  =  surf 
return 
endif 

if (x.gt .xl .and.x. It .x2) then 
c  to  calculate  the  area  of  the  grain. 

a  =  sqrt3  /  3.  *  (w  -  2 .  *  x  +  (p  +  2.  *  x)  /  2.)  **  2  -  pi  / 

1  6.  *  (w  -  2.  *  X  +  (p  +  2.  *  X)  /  2.  )  **  2 

theta  =  acos(dpc  /  (p  +  2.  *  x) ) 

al  =  theta  /  4.  *  (p  +  2.  *  x)  **  2  -  dpc  /  4 .  *  sqrt ( (p  +  2. 
1  *  x)  **  2  -  dpc  **  2) 

omega  =  acos(2.  *  dpc  /  (p  +  2.  *  x)  -  1.) 

a2  =  0.125  *  (p  +  2 .  *  x)  *  ( (p  +  2 .  *  x)  *  (omega  + 

1  sin (omega))  -  2.  *  dpc  *  sin (omega)) 

c  to  calculate  the  volume  of  the  grain. 

v2  =  6.  *  (gl  -  2.  *  X)  *  (2,  *  dpc  **  2  -  dpc  *  (p  +  2.  *  x) 
1  -  pi  /  4 .  *  (p  +  2 .  *  X)  **  2  +  2.  *  al  +  4 .  *  a2) 

c  to  calculate  the  surface  area  of  the  grain. 
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s2  =  2 .  *  v2  /  (gl  -  2.  *  x)  +  12.  *  (gl  -  2.  *  x)  *  (dpc  +  (p 
1  +  2.  *  x)  *  (pi  /  2.  -  omega  -  theta  -  sin (omega))) 

c  to  calculate  the  volume  and  the  surface  area, 
vf  =  v2 
sf  =  s2 

c  to  calculate  the  the  mass  fraction, 
frac  =  1  -  vf  /  VO 
surf  =  sf 

dsdx  =  -  surf  /  (x2  -  x) 
vol  =  vf 

dsdx  =  coefOcpr)  *  dsdxsl()cpr)  +  (1.  -  coef()cpr))  *  dsdx 
dsdxsl()cpr)  =  dsdx 

frac  =  coef  (Vtpr)  *  fracsl()cpr)  +  (1.  -  coef()cpr))  *  frac 
fracsl(kpr)  =  frac 

surf  =  coef()cpr)  *  surfsl()cpr)  +  (1.  -  coef()cpr))  *  surf 
surfsl(kpr)  =  surf 
return 
endif 

if (x . gt . x2) then 
dsdx  =  0 . 
surf  =  0. 
vol  =  0 . 

frac  =  fracsl(kpr)  *  coef (kpr)  +  1.  -  coef(kpr) 

fracsl(kpr)  =  frac 

if (frac .gt .. 9999)  frac  =  1. 

if (frac. gt. .9999)return 

dsdx  =  0. 

dsdx  =  dsdxsl(kpr)  *  coef (kpr) 
dsdxsl(kpr)  =  dsdx 
if (abs (dsdx) . It . 1 . ) dsdx  =  0. 
surf  =  surfsl(kpr)  *  coef (kpr) 
surfsl(kpr)  =  surf 
return 
endif 
stop 
c 

c  spherical  grain  calculations  start  here 
c 

210  if  (gd  .le.  2.*x)  goto  80 

vol  =  pi  /  6.  *  (gd  -  2.  *  x)  **  3 
surf  =  pi  *  (gd  -  2.  *  x)  **  2 
frac  =  ( (gd  -  2.  *  x)  /  gd)  **  3 
dsdx  =pi*4.  *  (2.  *x-  gd) 
return 
c 

220  format  (lx, ' unacceptable  granulation') 
end 

subroutine  jint (btdia, btlen, x, y, nchpts, chdist , chdiam, bint , bvol) 
dimension  bint (10) , chdist (6) , chdiam(6) 
pi=3. 141593 

areaa=pi*btdia*btdia/4 . 

distbp=x-btlen 

points=100 . 

step=y/points 

zz=0 . 

bvoll“0 . 

ichg=0 

do  1  j=l,10 
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bint ( j) =0 . 

1  continue 

if ( step . It . 1 . e- 10) then 

bint (7)=1./ (pi* (chdiam(l) **2/4 . ) ) **2 

return 

endif 

intsw=0 

do  2  j=2,nchpts 

if (y . gt . chdist { j) )  go  to  2 

nchp= j 

ichg=l 

holddt=chdist ( j) 
holddin=chdiam  ( j ) 

diarn- (y-chdist ( j-1) ) / (chdist ( j) -chdist ( j-1) ) 
chdiam( j) =chdiam( j-1) +diam* (chdiam( j) -chdiam( j-1) ) 
chdist ( j) =y 
go  to  3 

2  continue 
nchp=nchpts+l 
chdist (nchp) =y 
chdiam(nchp) =chdiam (nchpts) 

c  write (6, *) chdist (nchp-1) , chdist (nchp) 

c  910  format ( lx, 2ell . 4 ) 

3  continue 

areal=chdiam(l) **2  *  pi  /  4 . 

bint5o=0 . 0 

do  58  il=l, nchp-1 

npt=int ( (chdist (il+i) -chdist (il) ) /step) 
if (npt . le . 0 ) npt=l 

c  wr ite (6,  *) npt, chdist (il) , chdist (il+1) , step 

stepl= (chdist (il+1) -chdist (il) ) /npt 
c  write (6, 912 ) chdi am (nchp-1) , chdiam(nchp) 

c  912  f ormat ( lx, 2el 1 . 4 ) 

rl=chdiam (il) * . 5 
c  areal=pi*rl*rl 

do  57  i=l,npt 
zz=zz+stepl 

if (zz .gt . chdist (nchp) ) zz=chdist (nchp) 
diam= (zz-chdist (il) ) / (chdist (il+1 ) -chdist (il) ) 
diam=chdiam(il) +diam* (chdiam(il+l)-chdiam(il) ) 
r2=0 . 5*diam 
area2=pi*r2  *r2 

bvol2=bvoll+stepl*pi/3 .*(rl*rl+rl*r2+r2*r2) 
if (zz .gt .distbp) then 
if ( intsw . eq . 0 ) then 

diam= (distbp-chdist (ii) ) / (chdist (il+1) -chdict (i 1 ) ' 
diam=chdiam (il ) +diam* (chdiam (il+1 ) -chdiam (il) ) 
r2a=0 . 5*diam 
area2=pi*r2a*r2a 
stepla=stepl- (zz-distbp) 

bvol2=bvoll+stepla*pi/3 . * (rl*rl+rl*r2a+r2a*r2a) 
bintl=bint (1) +0 . 5*stepla* (bvoll/areal+bvol2/area2 ) 
bint (2) =bvol2*bvol2/area2/area2 

bint (3) =bint (3) +0 . 5*stepla* (bint (1) *areal+bintl*area2 ) 
bint (4)=bint (4) +. 5*stepla* (bvoll*bvoll/areal+bvol2*bvol2/area2) 
c  bint (5) =bint (5) + . 5*stepla* (1 . / area 1+1 . /area2 ) 

bint (6) =bvol2/area2/area2 
bint  (7)  =--l . /area2/area2 
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bxnt  (1) =bintl 

bvoll=bvol2 

areal=area2-areaa 

area2=pi*r2*r2 

area2=area2-areaa 

stepla=zz-distbp 

r l=r2a 

bvol2=bvoll+stepla*pi/3 . * (rl*rl+rl*r2+r2*r2) 
bvol2-bvol2-stepla*areaa 

bintlO=bint (10) +0 . 5*stepla* (bvoll/areal+bvol2/area2) 
bint5“bint (5)+0.5*stepla*(l. /areal  +1 . /area2) 
bint (2) =bvol2 *bvol2/area2/area2 
c  bint (3) =bint (3) +0 . 5*stepla* (bint (1) *areal+bintl*area2) 

bint (4) =bint (4 ) + . 5*stepla* (bvoll^bvoll/areal+bvol2*bvol2/area2) 
c  bint (5) =bint (5) + . 5*stepla* (1 . /area 1+1 , /area2) 

bint (6) =bvol2/area2/area2 
bint (7) =1 . /area 2 /area 2 

c  bint5a=bint5o+ . 5*stepla* (1 . /areal  +  l./area2) 

bint (8) =bint (8) +. 5*stepla* (areal*bint (5)  +  bint5*area2) 
bint ( 9) =bint ( 9) + . 5*stepla* (areal*bint (10) +area2*bint 10 ) 
bint ( 10 ) =bintl0 
bint ( 1 ) =bintl 
bint (5) =bint5 
c  bint5o=bint5a 

areal=area2 
bvoll=bvol2 
rl=r2 
intsw=l 
go  to  57 
else 

area2=area2-areaa 

bvol2=bvol2-stepl*areaa 

bint 10=bint ( 10) +0 . 5*stepl* (bvoll/areal+bvol2/area2) 
bint (2) =bvol2*bvol2/area2/area2 
c  bint (3) =bint (3) +0 . 5*stepl* (bint (1) *areal+bintl*area2 ) 

bint (4) =bint (4 ) +0 . 5*stepl* (bvoll*bvoll/areal+bvol2*bvol2/ area2 ) 
c  bint (5) =bint (5) + . 5*stepl* ( 1 . /areal +1 . /area2) 

bint ( 6) =bvol2/area2/area2 
bint (7) =1 . /area2/area2 

bint5=bint ( 5) + . 5*stepl* (1 . /areal  +  l./area2) 
bint (8) =bint (8) + . 5*stepl* (areal*bint (5)  +  bint5*area2) 
bint (9) =bint (9) +. 5*stepl* (areal *bint (10) +area2*bintl0 ) 
bint ( 10 ) =bint 10 
c  bint (1) =bintl 

bint (5) =bint5 
c  bint5o=bint5a 

areal^area2 
rl=r2 

bvoll=bvol2 
go  to  57 
endif 
endif 

bintl=bint (1) +0 . 5*stepl* (bvoll/areal+bvol2/area2) 
bint (2) =bvol2*bvol2/area2/area2 

bint (3) =bint (3) +0 . 5*stepl* (bint (1) *areal+bintl*area2) 
bint (4) =bint (4) +0 . 5*stepl* (bvoll*bvoll/areal+bvol2*bvol2/area2) 
c  bint (5) =bint (5) + . 5*stepl* ( 1 . /areal +1 . /area 2) 

bint (6) =bvol2/area2/area2 
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bint (7) =1 . /area2/area2 
bint ( 1 ) =bintl 
areal='area2 
r  l=r2 

bvcl l*bvol2 

57  rf'^cinue 

58  continue 
bvol=bvol2 

if (ichg.eq. 1) then 
chcliam(nchp)  =holdclm 
chdist (nchp) =holddt 
endif 

c  write(6,9l5) 

c  915  format (' 1' , lx, ' Leaving  Jint') 
return 
end 
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LIST  OF  SYMBOLS 


P 

A(z) 

A(za) 

A(za') 

ACza'") 

Aa 

Aba 

C 

8o 

m 

Mp 

P 

Pb 

Pbr 


P(za) 

t 

u 

V(x,t) 

V(za,t) 

V(zp,t) 


density 

area  at  position  z 
area  at  position  za 

area  just  before  the  aft  end  of  the  projectile 

area  just  after  the  aft  end  of  the  projectile 

cross-sectional  area  of  the  bore 

cross-sectional  area  of  the  boattail 

area  at  the  base  of  the  projectile 

total  charge  mass 

a  constant  to  reconcile  dimensions 

mass  flux 

mass  of  projectile 

pressure 

projectile  base  pressure 
breech  pressure 
mean  pressure 
resistive  pressure 

pressure  on  the  aft  end  of  the  projectile 

time 

velocity 

volume  as  a  function  of  distance  and  time 
volume,  at  za,  as  a  function  of  time 
volume,  at  zp,  as  a  function  of  time 
projectile  velocity 
derivative  with  respect  to  lime 
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ATTN:  SMCAR-FSA-T,  M.  Salsbury 

SMCAR-FSA-F.  LTC  R.  Riddle 
SMCAR-F,SC,  G.  Ferdinand 
SMCAR-FS.  T.  Cora 
SMCAR-FS-DH.  J.  Feneck 
.SMCAR-FSS-A,  R.  Kopmann 
B.  Machek 
L.  Finder 

SMCAR-FSN-N,  K.  Chung 
Picatinny  Arsenal,  NJ  OTSOb-.'ilXX) 

3  Director 

Benct  Weapons  Laboratories 
ATTN;  SMCAR-CCB-RA. 

G.P.  O'Hara 
G.A.  Pllegl 

SMCAR-CCB-S,  F.  Heiser 
Watervliet,  NY  12189-4050 

2  Commander 

U.,S.  Army  Research  OlTice 
ATTN;  Technical  Library 
D.  Mann 
P.O.  Box  12211 

Rc.search  Triangle  Park,  NC  27709-221 1 

I  Commander,  U.SACECOM 

R&D  Technical  Library 
ATTN:  ASQNC-ELC-LS-L-R. 

Mycr  Center 

Fort  Monmouth.  NJ  07703-5301 


1  Commandant 

IJ..S.  Army  Aviation  School 
ATTN;  Aviation  .Agency 
Fort  Rucker.  AL  36360 

1  Program  M;m;igcr 

U.S.  Tank-Automotive  Command 
ATTN:  AMCPM-ABMS.  T.  Dean 
Warren.  Ml  48092-2498 


Oreani/ation 


I  Project  Manager 

U.S.  Tank-Automotive  Command 
Fighting  Vehicle  Systems 
ATTN;  SFAE-ASM-BV 
Warren.  MI  48397-5000 

1  Project  Manager.  Abrams  Tank  System 

ATTN:  SFAE-ASM-AB 
Warren.  MI  48.397-.5000 

1  Director 

HQ.  TRAC  RPD 

ATTN:  ATCD-MA 

Fort  Monroe,  V.A  23651-5143 

I  Commander 

U.S.  Army  Bcivoir  Research  and 
Development  Center 
.ATTN:  STRBE-WC 
Fort  Belvoir.  VA  22060- .5(K)6 

1  Director 

U.S.  Army  TRAC-Ft.  l.ee 
ATTN:  ATRC-L.  Mr.  Cameroti 
Fort  Lee.  VA  2.3801-6140 

1  Commandant 

U.S.  Army  Command  :md  General 
SlalT  College 

Fort  Leavenworth.  KS  66027 

1  Commandant 

U.S.  Army  Special  Warlare  School 
A  I'TN:  Rev  and  Ting  Lit  Div 
Fort  Bragg.  NC  28307 

1  Commander 

Raillord  Army  .Ammunition  Plant 
ATTN:  SMCAR-Q.A/HI  LIB 
Radford.  VA  24141-0298 

1  Commander 

U.S.  Army  Foreign  Science  and 
Technology  Center 
ATTN:  A\LX.ST-,MC-3 
220  Seventh  Siieet,  NL 
Charlottesville.  VA  22901-5396 
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4 
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Ci)mmand;ini 

LI.S.  Army  Field  Arlillcry  Center  and 
Sch(H)l 

ATTN:  ATSF-CO-MW.  E.  Dublisky 
ATSF-CN.  P.  Gross 
Ft.  Sill.  OK  T.TSO.^-.SWK) 

Commandant 

U.S.  Army  Armor  .School 

ATTN:  ATZK-CD-MS.  M.  Falkovitch 

Armor  Ajiency 

Fort  Knox.  KY  4()l21-.S21.‘i 


Commander 

Naval  Sea  Systems  Command 
A1TN:  SEA  62R 
SEA  64 

Washington.  DC  20.^62-.Sl()l 
Commander 

Naval  Air  Systems  Command 
ATTN:  AIR-9.S4-Teeh  Library 
Washington.  DC  20-^60 

Commander 

Naval  Research  Laboratory 
.ATTN:  Technical  Library 
C  -Ic  4410. 

K.  Kaila.sanalc 
J.  Boris 
E.  Oran 

Was/iingiod.  tX'  20M5-5(H)i) 


3  Commander 

Naval  Surface  Warfare  Center 
ATTN:  CiKle  730 
Code  R-13. 

R.  Bernccker 
H.  Sandusky 

Silver  Spring.  MD  20903-5(KX) 

7  Commander 

Naval  Surface  Warfare  Center 
ATTN:  T.C.  Smith 
K  Rice 

S.  Mitchell 
S.  Peters 
J.  Consaga 
C.  Goi/mer 
Technical  Library 
Indian  Head.  MD  20f>4(i-.S(XK) 

4  Commander 

Naval  Surface  W;irlare  Center 
ATTN:  C(Hie  Ci.^0.  Guns  &  Mu  ni.ons  Div 
Code  Ci32.  Guns  Systems  Div 
Ctxle  G33.  T.  Doran 
Code  E23  Technical  Library 
Dahlgrcn.  VA  22448-.S(XK) 

5  Commander 

Naval  Air  Warfare  Center 
ATTN:  C(Hle  38S. 

C.F.  Price 

T.  Boggs 
Code  3«i).‘i. 

T.  Parr 
R.  Derr 

Information  Science  Division 
Chma  Lake.  CA  y3.‘).YV600l 


(.)fficc  of  Naval  Research 
ATTN:  Code  473.  R.S.  Miller 
S(K)  N.  Quincy  Street 
Arlington.  V A  22217-9999 


Commanding  Officer 
Naval  Underwater  Systems  Center 
ATTN:  Crxie  .SB331.  Technical  Library 
Newport.  Rl  02840 


Office  of  Naval  Technology 
ATTN:  ONT-2 1 3.  D.  Siegel 
8{K)  N.  Quincy  St. 

Arlington.  VA  222I7-.S(XX) 


I  AFOSR/NA 

ATTN:  J.  Tishkoff 

Bolling  AFB.  D  C.  20332-M48 

I  OLAC  PL/TSTL 

ATTN;  D.  Shiplctt 

Edwards  AFB.  CA  93.S23-.S(XX) 
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AL/LSCF 
ATTN:  J.  Levine 
L.  Quinn 
T.  Edwards 

Edwards  AFB,  CA  93523-5000 

WL/MNAA 

ATTN;  B.  Simpson 

Eglin  AFB.  FL  32.542-5434 

WL/MNME 

Energetic  Materials  Branch 
2.306  Perimeter  Rd. 

STE  9 

Eglin  AFB.  FL  32.542-.5910 

WL/MNSH 

ATTN:  R.  Drahe/uk 

Eglin  AFB.  FL  32.542-.5434 


1  Director 

Sandia  National  Laboratories 
Energetic  Materials  &  Fluid  Mechanics 
Department.  1512 
ATTN:  M.  Baer 
P.O.  Box  5m) 

Albuquerque.  NM  87185 

1  Director 

Sandia  National  Laboratories 
Combustion  Research  Facility 
ATTN.  R.  Carling 
Livermore.  CA  945.51-0469 

I  Director 

Sandia  National  Laboratories 
.MTN:  8741.  Cl.  A.  Beneditti 
P.O.  Box  969 

Livermore.  CA  94.5.51-0969 


NASA  Langley  Research  Center 
ATTN:  M.S.  408, 

W.  Scallion 
D.  Witcofski 
Hampton.  VA  2.3605 

Central  Intelligence  Agency 
Office  of  the  Central  References 
Dissemination  Branch 
R(X)m  GE-47,  HQS 
Washington.  DC  20502 

Central  Intelligence  Agency 
ATTN:  J.  Backofen 
NHB.  Room  .5N()1 
Washington,  DC  2050.5 


2  Direclor 

Lawrence  Livermore  National 
Laboratory 
ATTN:  L-.i55. 

A.  Buckingham 
M.  Finger 

P.O.  Box  808 

Livermore.  CA  94.55()-()622 

2  Director 

Los  .Mamos  Scientific  Lab 
ATLN:  T3/D.  Butler 

M.  Division/B.  Craig 
P.O.  Box  1663 
Los  Alamos.  NM  87544 


SDIO/TNI 

ATTN:  L.H.  Caveny 
Pentagon 

Washington.  DC  20.301-7100 

SDIO/DA 
ATTN:  E.  Gerry 
Pentagon 

Washingtoti,  DC  2 1.30 1-7  UK) 

HQ  DNA 
ATTN;  D.  Lewis 
A.  Fahey 

6801  Telegraph  Rd. 
Alexandria.  VA  22310-3.398 


2  Battelle 

.A 'Ll  N:  TACTEC  l.ibrary,  J.N.  Huggins 
V.  l.cvm 
.505  King  Avenue 
Columbus.  OH  43201-2693 

1  Battelle  PNL 

ATTN;  Mr.  Mark  Garnich 
P.O.  Box  999 
Richland.  WA  99352 

1  Institute  of  Gas  Technology 

AT  TN:  D,  Gidasfiow 
3424  S.  Stiite  Street 
Cliic;igo.  IL  60616-3896 
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1  Insiiiuic  lor  .Atlvaiiccd  Tcchiiolujiy 
ATTN:  T.M.  Kichne 

The  lliiivcrsiiy  of  Texas  of  Austin 
4010-2  W.  Riaker  Lane 
Austin.  TX  787.Sy-.-S.^29 

2  CPIA  -  JHLI 

A7TN:  H.  J.  Hoflntan 
T.  Christian 

10670  l.itlle  Patuxent  Parkway 
Suite  202 

ColiiinhKi.  ML)  2I(M4-.72(K) 

I  Briehani  Youiie  University 

Department  of  Chemieal  Engineering 
AriN:  .M.  Beekstead 
ProNO.  UT  84601 

1  .let  F’ropulsion  Lahoratory 

Calilornia  Institute  of  Teehnologv 
ATTN:  L.D.  .Strand,  MS  I2.S/224 
48(K)  Oak  Grove  Drive 
f'asadeiia,  CA  91100 

1  California  Insiuule  of  Technology 

204  Karman  Lah 
Main  Stop  .701-46 
,\7TN:  F.E.C.  Culick 
1201  E.  California  Street 
Pasadena.  CA  91109 


1  University  of  Maryland 

ATTN:  Dr.  J.D.  Anderson 
College  Park.  MD  20740 

I  University  of  Massrichusetts 

Depariinent  ol  Mechanieal  Engineering 
ATTN:  K.  Jakus 
Amherst.  MA  01(K)2-0014 

I  University  ol  Minnesota 

Department  of  Mechanical  Engineering 
ATTN:  E.  Fletcher 

Minnea(iolis.  MN  .S.S4 14-7.768 

7  Pennsylvania  State  Lniversity 

Depariiiiem  ol  Mechanieal  Engineering 
AITN:  Vang 

K.  Kuo 
C.  Merklc 

Umversity  Park.  PA  16802-7.701 

I  Rensselaer  Polylechnie  Institute 

rX'parliuciii  of  Maliieii.atics 
Troy.  NY  12181 

I  Stevens  Institute  of  Technology 

Davidson  Lahoratory 
AITN:  K,  McAlevy  111 
Castle  Poitit  Station 
Hoixiken.  NJ  07070-.s907 


7  Georgia  Institute  of  Technology 
SchiHil  of  .Aerospace  Engineering 
ATTN:  B.T.  Zim 
E.  Puce 
W.C.  Strahle 
.Atlanta.  G.A  70772 


1  Rutgers  University 

Departmeni  of  Mechanical  and 
Aerospace  Engineering 
AITN:  S.  Tcmkin 
Univeisity  Heights  Campus 
New  Bruns\Mck.  NJ  08907 


Massachusetts  institute  of  Technology 
Departmetit  of  Mechatitcal  Etigitieermg 
ATTN:  T.  Toong 
77  Massachusetts  .Avenue 
Cambridge.  MA  021.79-4707 

University  of  Illinois 
Department  of  Mechanical/Industry 
Engineering 
ATTN:  H.  Krier 
R.  Beddini 

144  MEB:  1206  N.  Green  St, 

Urbana.  IL  61801-2978 


I  University  of  Southern  California 

.Mechanical  Engineering  Department 
ATTN:  0HE2(K).  M.  Gerstein 
Los  Angeles.  CA  9(X)89-.7199 

1  University  ol  Utah 

Department  of  Chemical  Engineering 
ATTN:  A.  Baer 

■Salt  Lake  City.  UT  84112-1194 

I  Washington  Slate  University 

rX'partment  of  Mechanieal  Engineering 
.A  ITN;  C.T.  Crovvi- 
fhillman.  WA  99I6.7-.7201 
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1  AFELM.  The  Rand  Corporation 

ATTN:  Library  D 
1700  Main  Street 
Santa  Monica,  CA  90401-3297 

1  Arrow  Technology  Associates,  Inc. 
ATTN:  W.  Hathaway 

P.O.  Box  4218 

South  Burlington,  VT  05401-0042 

3  AAI  Corporation 
ATTN:  J.  Hebert 

J.  Frankie 
D.  Cleveland 
P.O.  Box  126 

Hum  Valley.  MD  21030-0126 

2  Alliaat  Tcchsy:  terns.  Inc. 

ATTN:  R.E.  lompkins 

J.  Kennedy 
7225  Northland  Dr. 

Brooklyn  Park.  MN  55428 

1  AVCO  Everett  Research  Laboratory 

ATTN:  D.  Stickler 
2385  Revere  Beach  Parkway 
Everclt.  MA  02149-5936 

1  General  Applied  Sciences  Lab 

ATTN:  J.  Erdos 
77  Raynor  Avc. 

Ronkonkaina,  NY  11779-6649 

1  General  Electric  Company 

Tactical  System  Department 

ATTN:  J.  Mand/y 
100  Plastics  Avc. 
r’illslicld.  MA  01 201 -.3698 

I  IITRI 

ATTN:  M.J.  Klein 
10  W.  35th  Street 
Chicago.  IL  60616-3799 

4  Hercules.  Inc. 

Radford  Army  Ammunition  Plant 
ATTN:  L.  Gizzi 

D.A.  Worrell 
W.J.  Worrell 
C.  Chandler 

Radford.  VA  24141-0299 


No.  of 

Copies  Ori’anization 

2  Hercules,  Inc. 

Allegheny  Ballistics  Laboratory 
ATTN:  William  B.  Walkup 
Thomas  F.  Farabaugh 
P.O.  Box  210 

Rocket  Center.  WV  26726 

I  Hercules,  Inc. 

Aerospace 

ATTN:  R.  Cartwright 
l(X)  Howard  Blvd. 

Kenville.  NJ  07847 

1  Hercules.  Inc. 

Hercules  Pla/a 
.\TTN:  B.M.  Rigg'  ‘man 
Wilmington.  DE  19X94 

1  MBR  Research  Inc. 

ATTN:  Dr.  Moshc  Ben-Reuven 
601  Ewing  St.,  Suite  C-22 
Princeton,  NJ  08.540 

1  01  in  Corporation 

Badger  Army  Ammunition  Plant 
ATTN:  F.E.’Wolf 
Barabw.  W1  5.3913 

3  Olin  Ordiumce 
ATTN:  E.J.  Kirschke 

A.F.  Cion/alez 
D.W.  Worthington 
P.O,  Box  222 

St,  Marks.  FL  32355-0222 

I  (7lin  Ordnance 

ATTN:  H.A.  McElroy 
lOIOI  9lh  Street.  North 
St.  Petersburg.  FL  33716 

1  Paul  Gough  Associates,  Inc, 

ATTN:  P.S.  Gough 
1048  South  St. 

Portsmouth.  NH  03801 -.5423 

1  Physics  International  Library 

ATTN:  H.  Wayne  Wampler 
P.O.  Box  .5010 

San  Leandro.  CA  94577-0599 
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2  PriiK'cion  Coinbusiioii  Research 

Luooralories,  Inc. 

ATTN:  N.  Mcr 

N.A.  Messina 
Pvinccion  Corporalc  Pla/a 
11  Deerpark  Dr.,  Bldg  IV.  Suite  119 
Monmouth  Junction.  NJ  08852 

}  RtK'kwell  International 

RiK'ketdyne  Division 
ATTN:  BAOK. 

J.  Flanagan 
J.  Gray 
R.B.  Edelinan 
66.2,2  Canoga  .Avenue 
Canoga  Park.  CA  91. 202-2702 

2  RtK'kwell  International  Science  Center 

ATTN:  Dr.  S.  Chakravarthy 
Dr.  S.  Palaniswamy 
1049  C'amino  Dos  Rios 
P.O.  Box  1085 
Thousand  Oaks,  CA  91260 

1  Science  Applications  International  Coip. 

ATTN:  M.  Palmer 
2109  Air  Park  Rd. 

.Alhuiiucrque,  NM  87106 

1  Southwest  Research  Institute 

ATTN:  J.P.  Ricgel 
6220  Culebra  Road 
P.O.  Drawer  28510 
San  Antonio.  TX  78228-0510 

1  Sverdrup  Technology.  Inc. 

ATTN:  Dr.  John  Deur 
2(X)1  Aerospace  Parkway 
Bnxik  Park.  OH  44142 

2  Thiokol  Corporation 
Elkton  Division 
ATTN:  R.  Wilier 

R.  Biddle 
Tech  Library 
P.O.  Box  241 
Elkton.  MD  21921-0241 

I  Veritay  Technology,  Inc. 

.^TTN:  E.  Fisher 
4845  Millcrsport  Hwy. 

East  Amherst.  NY  14.501-0,205 


No.  ol 

Copies  Oriiani/aiion 

1  Universal  Propulsion  Company 

ATTN:  H.J.  MeSpadden 
2.5401  North  Central  Ave. 

PhiK'inx.  AZ  85027-7827 

I  SRI  International 

Propulsion  Sciences  Division 
ATTN:  Tech  Library 
222  RavenwiKKl  Avenue 
Menlo  Park.  CA  9402.5-.2492 

.Alx.-rdceii  Provmc’  Ground 

I  Cdr,  USACSTA 

.ATTN:  STECS-PO/R.  Hendricksen 
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No.  of 

Copies  Ora;iiii/atioti 

I  Ernst-Mach-lnslitut 

ATTN:  Dr.  R.  Heiser 
Haupstras.se  18 
Weil  am  Rheim 
Gcnoaiiy 

1  Defence  Research  Agency.  Military 

Division 

ATTN:  C.  Woodley 
RARDE  Fort  Halstead 
Sevenoaks,  Kent.  TN14  7BP 
England 

1  Sch(H)l  of  Mechanical.  Materials,  and 

Civil  Engineering 
ATTN:  Dr.  Bryan  Lawton 
Royal  Military  College  of  Science 
Shrivenham,  Swindon.  Wiltshire.  SN6  SLA 
England 


Copies  Or»:iiii/:ilioi) 

2  Inslilul  Saint  Louis 

ATTN:  Dr.  Marc  Giraud 
Dr.  Gunther  Sheets 
Postlach  1260 
78.S8  Wcail  am  Rhein  I 
Germany 

1  Explosive  Ordnance  Division 

ATTN;  A.  W  ildcgger-Gaissmaier 
Dclence  Science  and  Technology 
Organisation 
P.O.  Box  I7.S0 

Salishury.  South  Australia  .SIOS 

1  Armaments  Division 

ATTN.  Dr.  J.  L;ivigne 

Defence  Research  E.stahlishment  Valcarlier 

24.Sy.  Pie  XI  Blvd..  North 

P.O.  Box  8800 

Couicclcite.  Oue'bec  GOA  IRO 
Canada 


III 


Intentionally  left  blank. 
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USER  EVALUATION  SHEET/CHANGE  OF  ADDRESS 


This  Laboratory  undertakes  a  continuing  effort  to  improve  the  quality  of  the  reports  it  publishes.  Your 
commenls/answers  to  the  items/questions  below  will  aid  us  in  our  efforts. 

1.  ARL  Report  Number  ARL-TR- 181 _ Date  of  Report  August  1993 _ 

2.  Date  Report  Received _ 

3.  Docs  this  report  satisfy  a  need?  (Comment  on  purpose,  related  project,  or  other  area  of  interest  for 

which  the  report  will  be  used.) _ 


4,  Specifically,  how  is  the  report  being  used?  (Information  source,  design  data,  procedure,  source  of 
ideas,  etc.) _ 


S.  Has  the  information  in  this  report  led  to  any  quantitative  savings  as  far  as  man-hours  or  dollars  saved, 
operating  costs  avoided,  or  efficiencies  achieved,  etc?  If  so,  please  elaborate. _ 


6.  General  Comments.  What  do  you  think  should  be  changed  to  improve  future  reports?  (Indicate 
changes  to  organization,  technical  content,  format,  etc.) _ 


Organization 


CURRENT  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 

7.  If  indicating  a  Change  of  Address  or  Address  Correction,  please  provide  the  Current  or  Correct  address 
above  and  the  Old  or  Incorrect  address  below. 


Organization 


OLD  Name 

ADDRESS  _ 

Street  or  P.O.  Box  No. 


City,  State,  Zip  Code 


(Remove  this  sheet,  fold  as  indicated,  tape  closed,  and  mail.) 
(DO  NOT  STAPLE) 


Department  of  the  Army 


OFRCIAl.  BUSINESS 


BUSINESS  REPLY  iU4IL 
HBi  ciAss  ffiiwr  >  oan,  Afs,  md 

Postage  *iH  be  psid  by  addressee 


Director 

U.S.  Army  Research  Laboratory 
ATTN;  AMSRL-OP-CI-B  (Tech  Lib) 
Aberdeen  Proving  Ground,  MD  21005-5066 


